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dynamics. Selected MATLAB programs are provided. 
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CHAPTER 1 
INTRODUCTION 


The concept of "Freedom" Space Station divulges a wide range of dynamics problems. 
The particular problem considered in this thesis will help examine how the mobile transporter 
motion affects the space station dynamics. In the current space station layout, the mobile 
transporter is connected to the main truss of the space station (see Figure 1). The mobile 
transporter moves along the length of the truss on a track carrying payload about the station. Since 
the sum of the mass of the mobile transporter and that of the payload it will carry is potentially 
comparable to the mass of the entire station, the inertial effects of the transporter should not be 
ignored. 

The current analysis is motivated by the Space Station-Mobile Transporter (SS-MT) 
system. A simplified model of a mass moving over a flexible guideway is used to resemble the 
more complicated SS-MT system. This simplified system may be solved with a continuous 
formulation, but obtaining a continuous formulation that will simulate the complicated system 
containing the station, transporter, and shuttle is not feasible. Therefore, while a continuous 
analysis provides insight, a more general discrete formulation is needed to address the large-scale 
problem at hand. The objective is to develop a discrete algorithm that can be tested against a 
continuous formulation for the Flexible Guideway/Moving Mass system. Once the validity of this 
discrete algorithm is proven, it may be extended to simulate the motion of the Space Station- Mobile 
Transporter system. 

In both the continuous and the discrete analyses, the guideway is modeled as a flexible 
beam. The moving mass is considered to be a rigid body that moves along the beam's length. The 
inertial effects due to the motion of the mass are included in the derivation. For the discrete 
system, the flexible beam is divided into finite beam-elements. Using lumped-parameter and 
energy-consistent methods, respectively, the beam's mass and stiffness matrices are obtained. 
Then, an invertible operator is used to map the continuous spatial representation of the moving 
mass onto the discrete representation of the flexible structure. 

A continuous and discrete formulation is used to examine the mass connected to the flexible 
guideway at only one point. This model does not correspond to the actual Space Station-Mobile 
Transporter system but it serves as an appropriate test article. A continuous formulation is 
developed to check the results of the discrete formulation, and numerical results are presented for a 
simply-supported and a free-free flexible beam. Once the validity of the discrete methodology is 


1 



Figure 1. 



; Station Assembly. 
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established, parametric studies are performed for both the simply-supported and the free-free 
beam. 


The final step in the analysis is to connect the moving mass to the flexible beam at two 
points. This model was chosen to represent the train/track aspect of the SS-MT system. Since a 
discrete methodology was validated, a continuous formulation of this new system is not necessary. 
To stay close to the SS-MT system, only the free-free beam case is examined here. Studies that 
examine how the spacing of the two contact points affects the dynamics of the entire system are 
presented. 

There are four chapters and four Appendices following this introduction. Chapter 2 
examines the relevant work that has been accomplished in this area. Most previous work 
concentrate on the motion of a truck travelling over a flexible bridge or on a train travelling over a 
track. The SS-MT dynamics are similar to the flexible bridge problem except for the rigid body 
motion that the inertially-free space station undergoes. Due to the increase in computer power in 
recent years, many dynamic simulation programs were developed. Chapter 2 examines how well 
these existing dynamic codes can handle the specific problem at hand. 

Chapter 3 describes the theoretical analysis: Section 3.1 focuses on the continuous 
formulation and Section 3.2 develops the discrete formulation. The simply-supported and the free- 
free beam are examined here as special cases. Section 3.3 discusses the discrete, inertially free, 
multipoint-of-contact system. Each of the three sections in Chapter 3 starts with a mathematical 
model of the system. Then the equation of motion is developed using a Newtonian method. The 
equations are made nondimensional and placed into a reduced set of modal coordinates. The final 
equations are then cast in state-space form, which is well suited for numerical computation. 

Chapter 4 displays the numerical results and is divided into two sections. The first section, 
Section 4.1, outlines how the simulations are organized and discusses the goals that the results are 
trying to obtain and the parameters used to develop the simulations. The second section, Section 
4.2, displays and discusses each simulation that is run. In this section, the simply- supported beam 
is examined first to compare the accuracy of the discrete and continuous formulation. Parametric 
studies are then performed in order to examine how a change in the system nondimensional 
parameters changes the system dynamics. Next, the free-free case is examined. Once again, the 
discrete and continuous formulations are compared to assess the accuracy of the discrete 
formulation. Parametric studies are again performed when the accuracy of the discrete case has 
been proved. Finally the two-point-of-contact case is displayed. 
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Chapter 5 offers conclusions and suggestions for further work; in particular, extending the 
discrete formulation presented here to model the Space Station-Mobile Transporter system. The 
concept of connecting the shuttle to the SS-MT system could also be considered, which would be a 
logical extension of the work presented in this thesis. 

Appendix A offers a Lagrangian formulation of the continuous free-free beam, which is 
used to check the continuous formulation of the same system derived earlier using Newton's 
equation of motion. Appendix B describes the numerical integration process that was used to 
obtain the simulation results shown in Chapter 4. The Runge Kutta integration scheme is 
described. Appendix C provides numerical values of the matrices used in evaluating the free-free 
beam and expands the spatial derivatives of the shape functions used in the discrete formulations. 
Appendix D describes the MATLAB programs that were written to simulate the Flexible 
Guideway/Moving Mass system. Some selected MATLAB programs are displayed. 
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CHAPTER 2 
PREVIOUS WORKS 


Reaching as far back as the early nineteenth century, scientists and theorists have been 
intrigued by the dynamic interaction that occurs when a load travels over a flexible structure. In the 
past, common models of this interactive system have been that of a track travelling over a bridge, a 
train travelling along a track, or a package moving on a conveyer belt. These systems are similar to 
the inertially free Flexible Structure/Moving Mass system that is used here to resemble the Space 
Station-Mobile Transporter (SS-MT) system. 

Unlike the other systems mentioned, the Space Station-Mobile Transporter system is 
inertially free, and it is geometrically complex. It is necessary to obtain a discrete representation of 
the Flexible Stracture/Moving Mass system that can be extended to the Space Station-Mobile 
Transporter system. A discrete representation is necessary in order to accommodate mass and 
stiffness matrices, rather than partial differential equations. References [1] and [2] were the first to 
address this important concept in a general way. This thesis is a detailed compilation of the 
analysis presented in those two papers, with the exception that Ref. [2] addresses the problem of a 
flexible structure moving on a flexible structure. 

This chapter introduces other work involving the dynamic interaction of the moving 
mass/flexible structure issue. Section 2.1 discusses some related previously released papers. 
Section 2.2 examines some existing multibody dynamics and finite-element codes in order to 
assess their ability to handle the proposed problem. 

2.1 PREVIOUSLY PUBLISHED PAPERS 

Three important features are needed to handle the Space Station-Mobile Transporter system 
and the dynamics of the Space Shuttle: 

(1) A discrete representation. 

(2) The ability to handle very large complex systems. 

(3) Computational efficiency. 

The previously published papers are generally not well suited to meet all three of these 
requirements. Most of the previous work assumes that a continuous model of the flexible structure 
is available (Refs. [3]-[7]). References [3]-[6] assume that the beam is simply- supported. For 
example, Galerkin's method (Ref. [3]), Inverse Laplace transform (Ref. [5]), and Fourier series 
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(Ref. [6]) are several of the analytical techniques discussed. The derived discrete methodology 
outlined in Section 3.2 uses the continuous solution presented in Ref. [5] for its comparison. 

Reference [7] presents a continuous formulation that may be applied to different boundary 
conditions. This continuous solution, however, is primarily beneficial when the inertial effects of 
the moving mass can be ignored, thus treating the travelling load as a moving force rather than as a 
moving mass. The repercussions of this assumption are disclosed in Chapter 4. 

Two papers (Refs. [8] and [9]) present discrete methodologies capable of solving the 
moving mass/flexible structure problem. Reference [8] presents a methodology that is applicable 
for time-varying forces. Reference [9] presents a formulation that may be implemented into a 
general finite-element code, such as MSC NASTRAN. This method is valid for any boundary 
condition. Lagrange multipliers are used to obtain a linear time-invariant formulation, which can 
then be solved using NASTRAN. However, as explained in Ref. [2], this formulation is generally 
applicable when modal reduction is not necessary. 

The papers discussed above (Refs. [3]-[9]) present methodologies that are well suited to 
study the motion of a heavy truck or train travelling over a flexible bridge. However, the 
methodology presented here is well suited for the complex inertially free Space Station-Mobile 
Transporter system. It is a discrete representation that is independent on the boundary conditions 
of the beam. It is well suited for model reduction so can be altered to include the Space Shuttle's 
dynamics. 

2.2 EXISTING CODES 

In the past twenty years, the availability and power of computers has grown exponentially. 
Coinciding with the boom in computer hardware was an increase in computer software capabilities. 
Today there are hundreds of software packages available to handle very diverse tasks. 

Among the codes that handle multibody dynamic interactions are: DISCOS (Ref. [10]), 
ADAMS (Ref. [11]), and DADS (Ref. [12]). These codes model the connection between two 
bodies as either a prismatic or a rotational joint. As was stated in Ref. [1], however, the motion of 
the Mobile Transporter is dependent on the flexible-body motion of the Space Station. Therefore, 
the SS-MT system cannot be modeled using these joints, without making far-reaching 
assumptions. 
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An alternative is to develop a formulation that could be implemented using an existing 
finite-element code, such as MSC NASTRAN (Ref. [13]) or ADINA (Ref. [14]). NASTRAN 
requires that the equations be cast in a time-invariant form. Ref. [9] presents an effective approach 
to using NASTRAN for the moving mass problem in the case of low-order structures. ADINA's 
strength lies in other areas. 

The discrete formulation presented in Section 3.2 is developed specifically with the 
application to large structural systems in mind. 
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CHAPTER 3 

THEORETICAL ANALYSIS 


Chapter 3 discusses the theoretical analysis: Section 3.1 presents the continuous 
formulation, Section 3.2 presents the discrete formulation with one point of contact, and Section 
3.3 presents a discrete formulation for the multipoint-of-contact system. 

Each section is set up in a similar manner. First, a mathematical model of the system is 
discussed. Using this model, equations of motion of the entire system are developed and placed 
into a nondimensional form in terms of the beam's modal coordinates. Then a state space matrix 
representation is formed so the formulation is suitable for numerical evaluation. 

A general formulation valid for any boundary condition is presented in each section. Two 
systems are examined in detail: (1) the simply-supported beam and (2) the free -free beam. 


3.1 CONTINUOUS FORMULATION 

In the following sections, the continuous formulation of the SS-MT system is developed. 
The results from this formulation are compared to those obtained by the discrete analysis presented 
in Section 3.2. 


Two different systems are examined here. The first system system model, shown in 
Figure 2, represents the inertially fixed space station with the transporter travelling along its 
length. This model represents many physical entities, the most popular example being a heavy 
truck travelling over a flexible bridge. The second system, shown in Figure 3, represents the 
inertially free space station with the transporter travelling along its length. 


A Xm 



Figure 2. Simply supported beam with a mass moving along its length. 
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Figure 3. Free free beam with a mass moving along its length. 

The dynamics analysis presented is a general methodology applicable to both systems. In 
Sections 3.1.5 and 3.1.6, the equations are made specific to the inertially fixed and the inertially 
free systems, respectively. 

3.1.1 Mathematical Model 

There are several ways to formulate the equations of motion for the systems depicted in 
Figures 2 and 3. For an inertially fixed system, an exact continuous solution using an inverse 
Laplace transform may be used (Ref. [5]). For a more general system, an exact continuous 
solution is not available. 

In this analysis, however, an assumed modes approach, which accurately gives the 
deformation while retaining only three modes, is used to solve the systems shown in the above 
figures. The assumed modes approach uses Galerkin's method to determine the beam deformation 
due to the motion of the mass. Galerkin's method assumes that the beam deflection can be 
approximated by a superposition of orthogonal mode shapes. The assumption is that, using 
orthogonal modes, the error of the approximation is vanishes as the number of modes retained is 
increased. 

The inertial effects of the moving mass are included in this analysis. If these effects are 
ignored, the mobile transporter is modeled as a moving force rather than as a moving mass. If the 
mass of the moving object is considered negligible compared to that of the flexible structure, the 
inertial effects can be ignored; however, if the mass of the object is comparable in size to the 
structure, these inertial effects should not be ignored. 
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Several cases of the two approaches for different moving object/flexible structure mass 
ratios were studied. The results are presented in Chapter 4. In the continuous and the discrete 
formulations, the flexible structure is modeled as a Bemouille-Euler beam. 

3.1.2 Equations of Motion 

The equations of motion for the flexible beam-moving mass system are derived in three 

steps: 

( 1 ) The partial differential equation describing the motion of the flexible structure is 
determined. 

(2) The moving mass equation is obtained. 

(3) A compatibility condition is invoked to obtain the equation of motion for the entire 
system. The compatibility condition is a direct consequence of Newton's third law. 

It states that for every force there is a reactive force equal in magnitude but opposite in 
direction. Thus, the force on the beam due to the moving mass is equal to the force 
created by the acceleration of the mass but opposite in direction. 

Flexible Structure Equation 


The fourth-order partial differential equation describing the motion for a Bemoulli-Euler 
beam is (Ref. [15]): 


p«M + £/^M =F ext (x,t) 

dx* (3.1) 

The displacement field, u(x,t), describes the motion of the mean chord of the deformed beam with 
respect to an inertial reference frame. This field contains any rigid body translation and rotation as 
well as flexible motion of the beam. Since the structure is modeled as a Bemouille-Euler beam, the 
mean chord of the beam is assumed to undergo pure translation as a result of the deformation. F ext 
is the vector sum of the external forces applied to the beam. The material properties of the beam 
are considered to be homogeneous. 



Moving Mass Equation 


The equation describing the motion of the moving mass is considered. This equation is 
derived using Newton's second law of mechanics, which states that the force exerted on the mass 
is proportional to the absolute acceleration of the mass. 

= F m (x m {xA‘) 

{ dt 2 labs (3.2) 


where F m (x m (x,t),t) is the force exerted on the moving mass. 


The absolute acceleration of the beam is the acceleration of the moving mass displacement 
field, u m (x m (x,t),t), with respect to an inertial reference frame. When differentiating the moving 
mass’s displacement field, it is important to imagine a reference frame embedded in the moving 
mass. Since the mass is moving, the reference frame is also moving with respect to the inertial 
frame and adds its own terms to the acceleration of the mass. In this analysis, the mass is moving 
at a constant velocity v m . Therefore, the absolute acceleration of the moving mass is 


d 2 M w j 

„ dt 2 Jabs 


— + 2 v m 


d 2 u m . 2 d 2 u m 
dtdx m m dxl 


(3.3) 


where the functional dependence of the variables are omitted for brevity. 

As it stands, Eq. (3.3) is written as a function of the moving mass displacement field, u m . 
To easily formulate the system equation, it would be helpful to express this equation in terms of the 
beam’s displacement u. Since the mass is fixed to the beam, it cannot slip and, at any time, the 
beam and the moving mass have the same displacement in terms of F m- This relationship can be 
expressed mathematically by using the dirac delta function, which is a continuous function that 
depicts the value of a function at discrete times. Using this, the displacement of the mass is 
rewritten as a function of the displacement of the beam: 

U-m = u m fyx - x m ) (3.4) 

Since the dirac delta function is not a function of time, the spatial derivatives in Eq. (3.3) 


are rewritten as 


(3.5) 



8 2 m 

dx 2 


S(x - x m ) 


d 2 u m _ d 2 u 
dx m dt dx dt 


- x m ) 


(3.6) 


Equations (3.5) and (3.6) are substituted into the absolute acceleration expression for the moving 
mass. The force equation for the moving mass then can be written in a form more compatible to 
the flexible beam equation: 


m m 


u + 2v m 


d 2 u 2 d 2 u 
dite +Vm dx 2 


8{x - x m ) = F m 


(3.7) 


Compatibility Equation 


Newton's third law of dynamics states that for every force there is a reactive force equal in 
magnitude but opposite in direction. Using this law, the compatibility equation between the 
moving mass and the flexible beam is determined. The force exerted on the beam due to the 
moving mass is equal in magnitude but opposite in direction to the force defined by Eq. (3.7). The 
total force exerted on the beam is the sum of the force due to the moving mass plus any other 
external forces acting on the beam: 


F ext ~ fext ~ F n 


(3.8) 


where fext is any arbitrary external force. The external force applied to the beam varies depending 
on which environment is being simulated. The actual value of fext f° r the two systems studied is 
shown in Sections 3.1.5 and 3.1.6. Using Eq. (3.7) the total force applied to the beam is 


F ext ~f ext ~ m n 


M + 2 v„ 


0 2 U 

dtdx 


+ v 2 


3 2 u 
dx 2 


S(x-x m ) 


(3.9) 


System Equation 

The equation of motion for the entire system is obtained by using the force expression 
found in Eq. (3.9) and substituting it into the beam equation (Eq. (3.1)). For clarity, any terms 
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involving the displacement of the beam are shown on the left-hand side of the equation even 
though they appear due to the force exerted by the moving mass. A fourth-order partial differential 
equation describing the total displacement of the beam as a mass moves at a constant velocity along 
its length is 


d*u 

pu + El — - + m m 

a * 4 


u + 2v„ 


a 2 u . 0 a 2 u 


dtdx 


+ vl- 


dx 2 


S(x - X m ) —f ext 


(3.10) 


The displacement field can contain translations and rotations. The boundary condition of 
the beam does not change the form of this general equation. 


Modal Solution 


Equation 3.10 must be rewritten into a form more suitable for numerical computation. For 
certain systems, this equation can be solved directly using an inverse Laplace transform. A more 
general approach is described here. A linear superposition of orthogonal modes is used to 
represent the beam's vibration. The new modal representation contains a mode shape, <f>i, which 
is only dependent on space, and a modal coordinate, 1 7 /, which is only a function of time. 

N 

u{x,t) = L X <Pi(x) ^ 

where N is the number of modes. As the number of modes approaches infinity, the modal 
approximation approaches the exact displacement of the beam u(x,t). The actual modes used 
depend on the system being analyzed. The mode shapes used for the simply- supported and free- 
free beam systems are shown in Sections 3.1.5 and 3.1.6, respectively. 

It is not feasible to use an infinite amount of modes to model the displacement. Different 
variations of the assumed modes approach that drive the error of the approximation to zero for a 
small amount of modes. have been developed. One of these methods is known as Galerkin s 
method, which uses the orthogonal property of the modes to drive the error to zero. 

Equation (3.1 1) is used to substitute the modal coordinate r\i(t) for the natural coordinate 
u(x,t)- this substitution is used in Eq. (3.10). This new equation can be integrated since the mode 
shapes (pi are not a function of time. Before the integration takes place, the new equation is 
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premultiplied by the transpose of the mode shapes <pj(t). Since the modes are orthogonal, this step 
reduces the error of the approximation. The resulting integral equation is 


r L 


N 


m m 


El fa rii + p fa rji j dx + 


fa 2 ifa + 2 v m fa r\i +v OT <Pi rii\@ Xm -faextj 




i = 1 


( 3 . 12 ) 


where 


frjextj ~ 



fa f ext 


(3.13) 


The roman numeral superscript indicates a spatial derivative of the appropriate order. A special 
relationship involving the dirac delta function was used in obtaining the above equation 


f(x) 8{x - x m ) dx =/(*«) 


(3.14) 


Nondimensional Integrals 

The mode shapes are only a function of space; therefore, the integral in Eq. (3.12) can be 
determined either analytically or numerically. To condense the equation into a more readable form, 
two nondimensional integrals are defined (Eqs. (3.15), (3.16)). The values of 1 1 and 12 for the 
simply-supported beam are found in Section 3.1.5. The values of the free-free beam are listed in 
Appendix C. 

h 0V)= f f fa fa dx 

Mo (3.15) 

I 2 ( ij ) = O f <pr fa dx 

Jo 


(3.16) 



Equation (3.12) is rewritten using the definitions in Eqs. (3.15) and (3.16). 


m m 


2 (p L 1 1 ( ij ) iii + jj 1 2 (ij>m }* 

(pj 2 (ft iii + 2 V m (pi Tji +V OT (pi T]ij@ x m = frjextj 




i= 1 


(3.17) 


This equation represents n first-order differential equations describing the total modal displacement 
of the system. Next, Eq. (3.17) must be put into a nondimensional form. 

3.1.3 Nondimensional Equations of Motion 

A general procedure for placing the equations of motion into a nondimensional form is 
presented. First, the reference parameters for each variable in the equation must be defined. Then, 
each variable is divided by the appropriate reference parameter. The reference parameters for the 
mass, time, and length variables are: 


Mass — » pL 

(3.18) 

Time — » -f- 

v m 

(3.19) 

Displacement — > L 

(3.20) 


The mass terms are made nondimensional by the beam's mass. The time reference parameter is the 
time required for the moving mass to travel the beam's length. The beam’s length is used as the 
reference parameter for displacement 


For clarity, nondimensional parameters are defined below and appear after the terms are 
divided by their respective reference parameters. 


Qi 


= (Oit r 


(»i L 

Vm 


(3.21) 


Mm 
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( 3 . 22 ) 


X 


(3.23) 


_ El t i = EL 

pL 4 ' pvjh L 2 

T = t - t\ m 

t r L (3.24) 

fi]extj 2 frjextj 

} = ^ 

P L P v « (3.25) 


where is the reference time defined above. The parameters represent the frequency, mass, 

* m 

stiffness, time, and external force of the system, respectively. In the following equations, an 
over-script o is used to represent a derivative with respect to the nondimen sional time parameter x, 
i.e.. 



Vm dt 


(3.26) 


Using the nondimensional procedure outlined previously, Eq. (3.17) is made 
nondimensional: 


X (h(iJ)°li + Xh(iJ)Vi\ + 

i= N l j = 1,2, ... N 

+ rn m <pj ^ TI i + 2 (p t Tji + <p 7]i)@ Xm = Hfj 

i-l (3.27) 


Equation (3.25) describes the displacement of the modal coordinate due to the motion of a moving 
mass and an externally applied force. Next, this equation must be rewritten in a form suitable for 
numerical simulation. 


3.1.4 State Space Representation 

The system depicted in Eq. (3.27) is described using a state-determined mathematical 
model. In this type of model, the system is described by a set of ordinary differential equations in 
terms of state variables (Ref. [16]). The future of all the variables associated with the system is 
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predicted from the previous time history of the state variables. The only information needed about 
the system is the initial condition of the state variables and the equations defining the future time 

history of these variables. 

To obtain a state space representation, an nth-order differential equation must be 
transformed into n first-order differential equations. Equation (3.27) is already in the required 
form. Next, an arbitrary set of state variables are chosen. For this model, the modal 
displacements and their associated velocities are chosen as the state variables: the displacements are 
chosen since they are the desired output, and the velocities are chosen so the matrices take on a 
familiar form. The two sets of state variables are combined into one state vector 

X (3.28) 

In a state-determined formulation, the time derivative of the state vector is a function of the 
state variables 


°x=F(x) (3.29) 

Eauations (3 28) and (3.29) show that the acceleration at any point can be expressed as a function 
X velolty and displacement at that point. Once the accelerations are known, the velocities and 
displacements are obtained by numerically integrating the system equation of motion forward m 

time. 


Matrix Representation 

For easy evaluation, Eq. (3.27) is placed into a matrix representation describing the states 
of the system. The equation is first rewritten so that it follows the standard matrix equation 

describing a dynamical system 

M 7] + C rj + K ri = F (3.30) 

where M, C, K, and F represent the mass, damping, stiffness, and force matrices, respectively. r\ 
is the vector containing the modal displacements 

T) = [i7i---nw] (3.31) 
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The definition of the state vector is used to place Eq. (3.30) into the form of Eq. (3.29). The final 
result is an equation that can be integrated to obtain the modal displacements and the modal 

velocities: 


0 

E 

X + 

0 

. -M l K 

-M a C . 


-M X F- 


(3.32) 


where E is the identity matrix. 

Additional terms representing structural damping are added to the damping matrix. It is 
easier to represent structural damping in modal coordinates rather than natural coordinates; 
therefore, the additional terms are already in modal form. Any off-diagonal modal terms are 
assumed to be negligible so the only extra terms appear on the diagonal. From Ref. [15], modal 

damping takes the form 

[Coh = 2 cVa Mm) I 2 (U) (3.33) 

where X2j is the nondimensional frequency of the beam. The beam's frequency depends on its 
boundary condition. The natural frequencies for the two systems examined are shown in Sections 
3.1.5 and 3.1.6. 

The mass, stiffness, and damping matrices all contain a constant and a time-varying 
component. The form of the force matrix depends on the type of external force applied to the 
system. The constant matrix is diagonal and represents the dynamics of the flexible beam without 
any moving mass. The time-varying matrix is fully populated and comes directly from the inertial 
effects of the moving mass. The combination of both matrices forms the total mass, stiffness, and 
damping matrices that are fully populated and time-varying. The constant matrices have the 
subscript o and the time- varying matrices have the subscript var . The total matrices are the sum of 

the two: 


M — Mq + A/ var 


(3.34) 


C = 


C 0 + c 


var 
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(3.35) 


(3.36) 


K = K 0 + K var 

First, the constant mass and stiffness matrices are determined. The constant modal 
damping matrix was defined by Eq. (3.33). Even though the actual values of the integrals have not 
been shown, the mode shapes used are orthogonal, ensuring the corresponding matrices will be 
diagonal. 


[m„] m = /,(/,() 

(3.37) 

My = A/ 2 (/,i) 

(3.38) 


Next, the time-varying matrices are shown. Note that in the following matrices the mode 
shapes and their derivatives are evaluated at the position of the moving mass, x m . Since x m 
depends on time, the values of the matrices also vary with time. 


[M va r]iJ — Hm $i 0/ 

(3.39) 

[C V ar]iJ ~ 2 /i m fa (pj 

(3.40) 

[Kvar]iJ = Mm 0/ $j 

(3.41) 


As stated previously, the vector containing the external forces may or may not be time 
varying, depending on the actual value of the external force applied. In symbolic terms, the force 
vector is 


F = 


L Mfy 


(3.42) 


The matrices shown in the Eqs. (3.37), (3.38), (3.39), (3.40), (3.41), and (3.42) are used 
in Eq. (3.32). The resulting expression is then numerically integrated to obtain the modal 
displacements and velocities at every point in time. Using Eq. (3.1 1) the modal displacements are 
transformed into the desired natural displacements. 
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To reiterate, the analysis presented so far has been for the general flexible beam/movmg 
mass system. Two systems are examined in detail using numerical methods: the inertially fixed 
system, which is modeled using a simply supported beam, and inertially free system, which is 
modeled using a free-free beam. Both systems are explained in greater detail in Sections 3.1.5 and 

3.1.6, respectively. 

3.1.5 Applications to a Simply-Supported Beam 

The simply supported system shown in Figure 2 represents many different physical 
systems. The most common physical system associated with this model is a truck traveling over a 
flexible bridge. The simply- supported beam model is used in this analysis to check the discrete 
methodology. One advantage of using this model is the availability of previously published 
results. Another advantage of simulating this system is the simplicity of the mode shapes. The 
nondimensional integrals can be calculated by hand; therefore, the numerical code used to simulate 
the system is easily checked when the simply-supported modes are used. 

The modes of a simply- supported beam are (Ref. [17]) 

<t>i{x) = sin l -&£ (3.43) 

These modes are orthogonal but not orthonormal. The corresponding nondimensional frequencies 
of the beam are 


A = = [infil (3.44) 

which are used in Eq. (3.33) to determine the constant component of the nondimensional damping 
matrix. 


Once the mode shapes are known, the values of the nondimensional frequencies are 
determined. For the mode shapes shown in Eq. (3.43), the integrals are determined analytically: 

7 'M = 2 (3.45) 


,. ., (i nf 
h (m) = '-y- 


(3.46) 
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The mode shapes shown in Eq. (3.43) are also used to develop the matrices given by Eqs. 
(3.33M3.41). Since the mode shapes for this system are simple sine waves, the matrices can 
easily be put into their symbolic form: 


M = 


2 + 

H m sin nr sin nr 


sym 


H m sin nt sin Nnr 


1 + 


fi m sin Nmtsin Nnrj 


(3.47) 


Cn 

*2 ••• fi m Nsinpr cosNnr 

+p m sinnt cosnr 


C =2 n 


fi m sinNnr cosnr 


£N 2 n 

+/x m /Vsrn Nnr cos Nnr 


(3.48) 


i- jt*X - 
2 

n 2 sinnr sinnr 


K = 


-n 2 sin Nnr sin nr 


-N 2 n 2 sin nrsinNnr 


t(N nfX- 

N 2 n 2 sinN nrsinNnr _ 


(3.49) 


It is easy to see the constant and time-varying components of these matrices. It is also 
apparent that only the total mass matrix and the constant components of the damping and stiffness 
matrices are symmetric. For more complex systems it is harder to write these matrices in their 
symbolic form. 
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to the beam is the gravitational force of the moving mass, .fat = -m m g Sfx ■ x m ). Since the mass 
is moving, this force varies with time. Using Eqs. (3.13) and (3.14), the force vector used for this 

simulation is 


F = 


fig sin Jt x 
fig sin Njct 


(3.50) 


where fig is a nondimensional parameter for the gravitational force applied to the beam: 

m m g .2 _ m m g 

^ 8 pL 2 P^m (3.51) 

The results of this simulation are presented in Chapter 4. 


3.1.6 Applications to a Free-Free Beam 

The free-free system, shown in Figure 3, may represent a crude model of the space station- 
mobile transporter system as it orbits around earth. The transporter is connected to the space 
station at one point. Therefore, this model depicts the transporter as a wheel travelling over the 
miss, rather than as a train travelling on a track. The train/track aspect of the mobile transporter is 

examined in Section 3.3. 

The model used to describe the menially free system is a free-free flexible beam with the 
rigid transporter travelling along its length. The modes of a free-free beam are (Ref. [17]) 

<t> 1=1 

02 =5c- 1/2 (3.52) 

(pi = cos pi x + cosh PiX -Oi (sin Pi x + sink Pi l) 3<i<N 


where 



(3.53) 
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(3.53) 



The first mode corresponds to rigid body translation. The second mode corresponds to 
rigid body rotation. The next N+2 modes are the flexible modes of the free-free beam. All the 
mode shapes are nondimensional, orthogonal, and orthonormal. The values for jBji and c; for the 
first three flexible modes are located in Appendix C. 

The corresponding nondimensional frequencies of the free-free beam are 

Qi = A 2 VX (3 .54) 

which, for this system, are the frequencies used when forming the constant modal damping matrix. 

For the free-free beam system, the nondimensional integrals are not determined 
analytically. Instead, // and 12 are determined by numerical integration. A fourth-order Runge 
Kutta integration scheme (detailed in Appendix B) is used, which is the same integration scheme 
used to integrate the equations of motion. 

The constant matrices are formed using the nondimensional integrals and the 
nondimensional frequencies. The time-varying matrices are formed using the mode shapes given 
in Eq. (3.52) at the appropriate value of T. There is nothing gained by writing out the specific 
matrices in their symbolic form for this system. The constant mass, stiffness, and damping 
matrices for the first three flexible modes are available in Appendix C. 

Since this system is designed to model an inertially free system, there is no gravitational 
field present. An initial vibration or an external force is needed to excite the system. For this 
simulation, an initial vibration was used rather than an external force. When the SS-MT system is 
attached to the shuttle system it is possible for the first mode to be excited due to the attitude control 
system of the shuttle. To create an initial excitation the left and right tip deformations were set 
equal to .02L with contributions from the first mode only. The moving mass was then released 
onto the beam as it was vibrating. The results are presented in Chapter 4. 
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3.2 DISCRETE FORMULATION 


Section 3.2 develops the discrete formulation of the SS-MT system that was analyzed in 
Section 3.1. The results obtained from this derivation are compared to the results obtained by the 
continuous formulation derived in Section 3.1. As before, two different systems are examined. 
The first system model, shown in Figure 2, represents a moving mass traveling along an mertially 
fixed structure. The second system, shown in Figure 3, represents the mass moving over an 
inertially free structure such as the space station. 

A general methodology is presented for analyzing any system. In Sections 3.2.5 and 3.2.6 
the methodology is made specific for the two systems described above. 

3.2.1 Mathematical Model 

In this formulation, the continuous systems of the previous sections will be placed into a 
discrete representation. As stated previously, a Bemouille-Euler beam is used to model the flexible 
structure. Discrete mass and stiffness matrices are determined for the flexible beam. The 
deflection of the beam and all the external forces applied to the beam are made discrete by 
introducing an invertible operator that distributes the effects of the moving mass over the 
appropriate discrete elements. 

First, the discrete equation of motion for the flexible structure is developed by creating 
discrete mass and stiffness matrices using either finite-element or lumped-parameter methods. 
These matrices represent the physical properties of the flexible structure. The goal is to discretize 
the load exerted on the beam due to the moving mass so it can be used with the already existing 
property matrices. To achieve this, a vector is formed that distributes the continuous forces along 
the beam's discrete points. A vector is created using two different finite-element shape functions: 
linear and cubic, which are compared in Chapter 4. Using the equivalent forces, the discrete 
equation is formed. To coincide with the continuous formulation the discrete equation is formed m 
terms of modal coordinates. This equation is then made nondimensional and placed into state 
space domain. The results for the simply-supported and the free-free beam are shown in Chapter 

4 . 
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3.2.2 Equations of Motion 


The discrete matrices that represent the mass and stiffness of the beam are determined and 
are used to write the general matrix equation of motion for a beam. This equation is the same one 
shown in Eq. (3.30) but is rewritten here for convenience 


Mq+Cq+Kq=F ( 3 . 55 ) 

where M, C, K, and F represent the mass, damping, stiffness, and force matrices, respectively. 
When this equation was used in Section 3.1, the matrices used represented the physical properties 
of the modes used to describe the motion. In the above equation, the matrices represent the 
discrete properties of the different beam elements used to model the beam. Even though the two 
equations have the same form, they represent two different systems. 

First, the discrete mass and stiffness matrices are formed. A discrete matrix for the 
damping is not developed in this subsection; however, a modal damping matrix is introduced in 
Section 3.2.4. The discrete stiffness matrix is developed using finite-element (or energy- 
consistent) techniques. Two different discrete mass matrices are formed. One matrix, developed 
from finite-element techniques, is used when it is important to keep the rotational inertias of the 
beam elements. The other mass matrix, formed using a lumped parameter model, is used when 
only the translational degrees of freedom of the beam elements are required. 

Next, a vector is developed that weights the continuous force due to the mass over the 
discrete beam elements. This vector is also used to discretize the deflection due to the moving 
mass. By combining the discrete property matrices and the discrete forces due to the moving 
mass, the discrete equation of motion for the system is formed. 

Mass Matrix 

The mass matrix for the beam is derived using both a lumped-parameter analysis and the 
energy-consistent finite-element method. When the linear shape function is used, the rotational 
degrees of freedom are statically condensed out of the mass and stiffness matrices; therefore, a 
lumped-parameter model is easily used. When the cubic shape function is used, each element's 
rotational degrees of freedom are needed; therefore, the mass matrix will be developed using the 
finite-element method. 
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Lumped-Parameter Model. The lumped-parameter method is appropriate only when 
the beam's material properties are homogeneous. In this particular analysis, this requirement is 
met; therefore, the model is valid. First, the beam is broken up into n finite elements. Then the 
mass of each element is distributed between the two neighboring nodes. In the case of the lumped- 
parameter model, the mass contribution at each node is half the mass of each element. The mass of 
each element is 


m e = pl e 


( 3 . 56 ) 


where l e is the length of each element and is defined as 



( 3 . 57 ) 


Since the material properties are continuous throughout the beam, the total mass at each node is the 
sum of the contributions from the two neighboring elements. The total mass at each node is 


m i = 2 m **2 

= m. 


m e 


2 <i<n 


( 3 . 58 ) 


Since the first and last nodes only feel the effects of one finite element, the mass contribution at 
those nodes is half the mass contribution at the inner nodes. 


Each node has a corresponding translation and rotation. Since the rotational inertias of each 
beam element are so small, the rotational degree of freedom can be eliminated from the stiffness 
matrix by using static condensation. 

When the linear shape function is used to distribute the force, only the translations at each 
node are important. Therefore, the mass matrix should only contain the translational degrees of 
freedom, which is accomplished fairly easily in a lumped-parameter mass matrix. The lumped 
mass matrix is diagonal with every other row, starting with the first row, corresponding to the 
translational degrees. The other rows correspond to the rotational inertias of each node. Since this 
will not be included in the mass matrix, the rotational inertias have not been shown. The final 
translational mass matrix has n+1 degrees of freedom and is in the following form 
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n+1 


M t = 


•.0 0 
0 mi 0 

0 O ' ■. 
n+i 


(3.59) 


This matrix is constant and discrete. 

Finite-Element Model. The matrix shown in Eq. (3.59) is used with the linear shape 
function. However, when the cubic shape function is used it is necessary to have access to both 
the translational and rotational degrees of freedom. It is possible to simply add the rotational 
inertias of the elements into the lumped-parameter model shown above. Instead, however, an 
energy-consistent mass matrix is developed. The finite-element approach is used to show another 
way to obtain a discrete mass matrix and is also used for the stiffness matrix. Each element of the 

finite-element mass matrix is (Ref. [15]) 


( L 

mtj — I Si Sj dx 

Jo 


(3.60) 


where si is a finite-element trial function. To correctly model a beam element, Hermites cubics are 
chosen for the trial functions because they have a continuous spatial second derivative (Ref. [15]). 
A trial function is needed for the deflection and the slope at each end of the element. Therefore, 
four trial functions for each element are needed. The four cubics are shown below. 


*='-^r +2 (tr 

(3.61) 

*-£- J (Sf + (ST 

(3.62) 

*- 3 ®- 2 (£f 

(3.63) 

kT3 

+ 

i 

II 

2* 

(3.64) 
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To determine the mass matrix for one element, the above trial functions are substituted into Eq. 
(3.60). This expression is then integrated to obtain the elemental mass matrix. The mass matrix is 
partitoned into four different matrices 

m e l e [ mu mn 1 

420 L m2i m22 J (3.65) 


where me and le are the mass 


and length of each element 


mu = 


156 
22 l e 


mi2 = 


54 
13 l e 


mi = 


54 

-13 l e 


, respectively. 

22 l e ' 

4 l} . 

-13// 

- 31 }. 

13 l e ' 

- 31 } . 


The four matrices are 

(3.66) 

(3.67) 

(3.68) 


m2 = 


156 
-22 l e 


- 22 le 

41} . 


(3.69) 


Next, the elemental mass matrices are combined to form the final global mass matrix. At this point 
all the degrees of freedom, translational and rotational, are present. Since the inner nodes connect 
two consecutive elements, the elemental mass matrices overlap. Therefore, the final global mass 
matrix is 


m = 


m 

420 n 


mn mi2 0 0 

mi m n + m2 m i2 o 

0 

0 0 . mu + m 2 2 

0 0 0 m2i 


0 

0 

0 

m ]2 

m 2 


(3.70) 


Stiffness Matrix 

The stiffness matrix is developed using finite-element techniques. When the linear shape 
function is used, only the translations at each node are required. Therefore, the rotational degrees 
of freedom are statically condensed out. The global stiffness matrix is developed the same way as 
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the energy-consistent mass matrix. Using the finite-element method, the elemental stiffness matrix 
is determined by (Ref. [15]) 



(3.71) 


Once again the cubics shown in Eqs. (3.61)-(3.64) are substituted into Eq. (3.71). After 
integration, the elemental stiffness matrix is obtained and, like the mass matrix, is also partitoned 
into four different matrices: 


k _ (£7) e [ leu 

ll L *21 


k]2 
h 2 . 


(3.72) 


where (EI) e is the elemental bending stiffness. The four matrices are 


hi 


kn 


hi = 


h2 



12 

6l e 

6l e 

4l 2 e 

-12 

6l e 

-6l e 

2li 

-12 

-6l e 

6l e 

2li 

12 

-6l e 

-6l e 

4l 2 e 


(3.73) 


(3.74) 


(3.75) 


(3.76) 


Next, the element stiffness matrices are combined to form the global stiffness matrix. At this point 
all the degrees of freedom, translational and rotational, are present. Since the inner nodes connect 
two consecutive elements, the stiffness matrices overlap. Therefore, the final global stiffness 
matrix is 
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kji ki2 0 0 0 

(FJ\ ^21 kjl + k22 kj2 0 0 

k= ^ o \ ■*. 0 

l e 0 0 *\ ku + k22 kl2 

0 0 0 k 2 i k 2 2 J (3.77) 

When the linear shape function is used, Eq. (3.77) is altered to condense out the rotational 
degrees of freedom. This reduced stiffness matrix is used in conjunction with the mass matrix 
shown in Eq. (3.59). For the cubic shape function, the matrix, as it stands in Eq. (3.77), is used 
with the similar mass matrix shown in Eq. (3.70) that contains both the translational and rotational 

degrees of freedom. 

The global stiffness matrix of Eq. (3.76) can again be partitoned into four separate 
matrices, k w k tr , k n , and k w The subscripts indicate either translational or rotational degrees of 

freedom. The partitoned stiffness matrix is 

£ _ (jgk kit ktr 

ll k rt krr\ (3.78) 

As stated previously, the rotational degrees of freedom are eliminated when using the linear shape 
function. The rotational degrees of freedom are statically condensed out. This is achieved by 
using the static matrix equation in Eq. (3.79): 

kttktr { v I _ 0 

krlkrrl \ q j 0 J (3.79) 

where v is a generic translational coordinate and 6 is a generic rotational coordinate. Solving for 6 
in terms of the translation, v, the reduced stiffness matrix becomes 

Kt = k tt - ktr kJr krt (3.80) 

Eq. (3.79) is a square matrix that is constant and has (n+1) degrees of freedom. When the 
rotational degrees are not eliminated, the stiffness and the matching mass matrix has (2n + 2) 
degrees of freedom. 
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Beam Equation 


Using the mass and stiffness matrices defined above the beam equation of motion is 
determined as 


Mq + Kq=F t (3.81) 

where M and K are generic discrete matrices representing the appropriate mass and stiffness 
matrices, depending on which case is being examined. The nodal displacements, contained in the 
q vector, represent the displacement at each node for the different elements. The nodal 
displacements, q, should not be confused with the modal displacements, 77 , discussed in Section 

3.1.2. The total discrete force vector. Ft, is a combination of any external forces applied to the 
beam and the inertial effects of the moving mass. This force vector is the discrete form of the 
vector F(. It correctly weights the effects of the moving mass onto the nodes of the beam. It is 

made discrete by using the discretization vector defined below. 

Discretization Vector 

In order to weight the effects of the continuous force between two discrete nodes, an 
invertible operator, called the discretization vector because it places the continuous forces into a 
discrete form suitable for Eq. (3.81) is developed. For any arbitrary time the discretization vector 
weights the effects of the moving mass. Since the mass is moving along the the beam, the vector 
must change with time to reflect this motion. Finite-element shape functions are used to distribute 
the forces. Two different shape functions are examined below. The first function, which is based 
on a linear interpolation, only looks at the translation at each node. The second function, which is 
of cubic order, takes into account the translation and rotation at both nodes. The difference 
between the two approaches is examined at length in Chapter 4. 

The two shape functions are developed in the same manner. A weighting function is used 
to locate the position of the moving mass with respect to the two appropriate nodes. The weighting 
function is defined as 


« = 


Xm ~ Xj 
Xi+1 - Xi 


(3.82) 


32 


The weighting function, l depends on the distance between two neighboring nodes, xi and x,+l. 
Note that Xi is defined as the nodal position that is either directly at or to the immediate left of t e 
moving mass. As soon as the point mass passes the Xi node location it is considered to be at die 
x i+ l position. The other variable, x m , has previously been defined as the location of the moving 

mass. 

The weighting function % in Eq. (3.82) is the discrete version of the continuous dirac delta 
function. The continuous function places the mass at a specific point, whereas $ weights the force 

between two neighboring nodes. 

To simplify the numerical evaluation of £, x[ is written in a suitable style for numerical 
evaluation. Two relationships are needed to accomplish this. First, as stated previously, the mass 
is considered to move at a constant speed; therefore, the position of the mass at an arbitrary time is 
always known. Next, it is assumed that the finite elements are of equal length. Using these facts, 
the distance between the two elements can be expressed as a function of the total beam length. 
These two relationships are shown symbolically as 


1! 

< 

3 

(3.83) 

Xi+l 

(3.84) 

Using the above two relationships, weighting function $ is rewritten as 


£ = ( v mf ' x i) jr 

(3.85) 


where j t/ is numerically calculated from 


Xi = int\^-\ m t)Ax ( 3 . 86 ) 

The int ( ) function truncates and retains only the integer portion of a real number argument. Note 
that if the i th node was defined as the node immediately to the right of the mass, then the position 
of xi would be rounded up rather than truncated down as shown in Eq. (3.86). 

Equations (3.85) and (3.86) are used to numerically evaluate the weighting function £. 
The first operator is linear with respect to the weighting function £ 
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Vl-{l-$V,+ SV M 


( 3 . 87 ) 


where the following vectors are defined as 


vf={o,o,---, o, l i,o, o of 


(l xf) 


vT +] = \ 0 , 0 0 , 1 ' , 0 , 0,--;0 


(mT 


\0 xf) 


( 3 . 88 ) 

( 3 . 89 ) 


where/is the degree of freedom for the particular system being analyzed. When the above vectors 
premultiply the vector of nodal displacements, velocities, or accelerations, they will locate the 
and (i+1 values respectively. In this method there are only translational degrees of freedom; 
therefore, only one value at each node is required. In the cubic shape function, however, it is 
necessary to capture two values at each node. 

When ^ is equal to zero, all the effects of the moving mass are placed at the i™ node. 
When £ is equal to one, all the effects are placed at the i^+7 node. For £ values between zero and 
one, the effects are appropriately weighted between the two nodes. 

The form of the linear interpolation function is easily determined without much 
computation. However, when the translations and the rotations at each node must be considered, 
the function's form is not easily seen. Therefore, the cubic shape function is developed in a more 
theoretical manner. 


Cubic Shape Function Definition 

As in the linear case, the cubic shape function is depicted as a function of £ but for the 
cubic function, the coefficients are defined in terms of £, and Also, in this case, there are 
two displacements at each node - translation and rotation. The values of each node can be thought 
of as the boundary conditions for the shape functions. When there are only two boundary 
conditions to be satisfied, a linear function will suffice. To satisfy four boundary conditions, 
however, a cubic function is required. 
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Since there are four values that need to be captured (two at each node), four V vectors are 
needed Vi Vi+ 1 Vi+2, and V/+J. The first and third vectors capture the translation at the i r and 
(i i + l)th nodes, respectively. The second and fourth vectors capture the rotation at the same 
respective nodes. Using these four vectors, the cubic shape function is determined using the 
standard finite-element method for determining shape functions, outlined below. 

In determining the cubic shape function, it is necessary to develop four trial functions that 
will multiply the four vectors described above (Ref. [18]): 

V 3 = Ti Vi + T i+ i Vi+i + T i+ 2 V i+2 + T i+3 V i+3 (3.90) 

Each trial function has the cubic form 

T, = a£ + + c,r + d£ 3 (3.91) 

The constants for each trial function are obtained by employing the four appropriate boundary 
conditions for each function (displayed in Table 1). 


Table 1. Boundary Conditions used in Determining Cubic Trial Functions. 


Trial Function 

S = o 

S = 1 

Ti 

ll 

II 

o 

Ti = 0 T# = 0 

Ti+1 

Ti+1 = 0 T&+l = l e 

Ti+1=0 Tfy+l = 0 

Ti+2 

Ti+2 = 0 T#+2=0 

Ti+2 = le Tfy+2 = 0 

Ti+3 

Ti+3 = 0 Tti+3 = 0 

Ti+3 = 0 Tti+3 = le 


In Table 1 the subscript £ indicates a derivative with respect to A similar table for the linear 
shape function could have been developed. However, in the linear shape function example it is 
trivial to develop the two trial functions. 


Using these conditions the four shape functions are determined. It is found that the 
appropriate trial functions are the Hermite’s cubics described in Eqs. (3.61)-(3.64). Substituting 
these trial functions into Eq. (3.90) leads to the final cubic shaping vector 
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V S = (i - 3 i 2 + 2 | J ) Vi + («- 2 | 2 + V lt , 
+ (j$ 2 -2{ 3 )v m + (-« J + « J )£Vm 


(3.92) 


For the general discrete derivation of the final equations, a generic discretization vector V is 
employed. This V vector represents either Vj or Vj, depending on which interpolation function is 

used. 


Load Modeling 

The force exerted on the beam is broken up into two components: the first encompasses 
any external force that is applied to the beam, and the second is the force exerted on the beam due 
to the inertial effects of the moving mass. The total force is the sum of both components. 

Ft = Fm + F m (3.93) 

Force Due to the Inertial Effects of the Moving Mass. The force due to the 
inertial effects of the moving mass is equal to the force created by the acceleration of the moving 
mass but it is opposite in direction. As seen in Section 3.1, the force is proportional to the absolute 
acceleration of the moving mass. This acceleration, however, is now given in terms of the discrete 
displacement field of the moving mass, q m . 

V dt 2 labs (3.94) 

where F m is the moving mass's discrete force vector. Equation (3.94) is the discrete counterpart 
of Eq. (3.2). In the following derivation, the functional dependence of the variables are omitted 
for brevity. 

The absolute acceleration of the discrete displacement field q m , must be determined. In 
order to do this, a relationship is needed between the beam's displacement field and the moving 
mass' displacement field. In essence, a discrete counterpart of Eq. (3.4) is needed, which is 
accomplished by using the discretization vectors defined above. For a general methodology, the 
generic shape function V is used. The relationship between qm and q is defined as 


q m = V T q 


(3.95) 
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Using Eq. (3.95), the absolute acceleration of the moving mass is written in terms of the 
beam's displacement. When determining the absolute derivatives, it is important to specify the 
variables of which V and q are a function. For simplicity it is assumed that V is a function of § 
only and £ is independently a function of time. The beam's displacement field q, is only a function 
of time. Using these conventions and the definition of an absolute acceleration determined in 
Section 3.1, Eq. (3.94) becomes 

F m = - m m V T q + 2 v\ \q + + V^($f ) < 7 ] (3.96) 

In this system, the mass is assumed to move at a constant velocity; therefore, £ is equal to zero. 
Using the definition of x shown in Eq. (3.85), its derivative with respect to time is 

c v m n 

1 = ^7 (3.97) 

Using Eq. (3.19), the above equation is rewritten as 

^ = t (3-98) 

where T was previously defined as the time required for the mass to move over the entire beam 
length. Written in this form, it becomes apparent that k represents a first-order discrete spatial 

derivative. For the duration of the general derivation, the spatial derivatives, V £, and v Q , are kept 
in their symbolic form. The actual values of both quantities for the linear and the cubic shape 
functions can be found in Appendix C. 

When using the linear shape function, the last term in Eq. (3.96) is zero. This term 
represents the force exerted on the beam when the mass moves over the beam's curvature. By 
examining the continuous case, it is apparent that this term adds a substantial force to the beam. 
Therefore, an impulse force is added to correctly model this force that results from the difference in 
slope of two neighboring elements (see Figure 4 and Section 3.2.6). 

When the higher-order, cubic shape function is used, the last term in Eq. (3.95) is not zero 
and the force from the beam's curvature appears without having to add an impulse force. 
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The inertial force due to the moving mass is now in a discrete form. The next challenge is 
to distribute the total force applied to the beam between appropriate nodes. Using Eqs. (3.93), 
(3.95,) and (3.98), the total force applied to the beam is 

Ft = F ex t - m m \v T q + 2 &V o\ (3 99 ) 

External Force. The external force applied to the beam is different for the two systems 
that are examined. For this derivation the external force is kept in its symbolic form. 

The total force is distributed between the appropriate nodes of the beam, which, in this 
analysis, is accomplished by using the same finite-element shape functions described in detail in 
the beginning of this section. Using these shape functions, the total discrete force is distributed to 
the appropriate nodes as 


F t = VF t 


(3.100) 


Using Eq. (3.99) this force becomes 


F t = VF ext -m m (w T q + 2(-f) VV^q + lff 

Equation (3.101) is the discrete form of the force shown in Eq. (3.9). 

System Equation of Motion 

The entire discrete equation of motion for the system is obtained by substituting Eq. 
(3.101) into Eq. (3.99) resulting in Eq. (3.102) below. As in the continuous case, the terms 
involving the beam's deflection q, are shown on the left-hand side of the equation, even though 
they appear due to the inertial effects of the moving mass. 

(M + m m V V T ) q + 2 m m * V v\q + (k + rn m V q^-VF^ 

Equation (3.102) is already in matrix form, unlike its continuous counterpart shown in Eq. (3.10), 
and is dimensional and in terms of the physical discrete beam coordinate q. Since the boundary 




(3.101) 
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conditions of the beam have not yet been specified, the above equation is valid for either the 
simply-supported or the free-free beam. 

3.2.3 Nondimensional Equations of Motion 

Unlike the equations in Section 3.1, Eq. (3.102) first is made nondimensional and then is 
transformed into the beam’s modal coordinates. For the discrete analysis, it is easier to perform 
modal reduction once the equation is in nondimensional form. Even though Eq. (3.102) is a 
matrix equation, the same general procedure is used to place the equation into a nondimensional 
form. A matrix is considered nondimensional if each of its elements are nondimensional. 
Therefore, each element is divided by the appropriate reference parameter, defined in Section 3.1 

When a common variable appears throughout an entire matrix, it can be extracted and 
placed in front of the matrix. This technique is used to define nondimensional mass and stiffness 
matrices. Referring to Eqs. (3.65) and (3.72), the new nondimensional matrices are 


M = J—M 

(3.103) 

pL 


(3.104) 


In a similar manner, the vectors containing the nodal accelerations, velocities, and 
displacements are made nondimensional by dividing each of their elements by the appropriate 
reference parameters. Again, some variables can be extracted and placed in front of the vectors. 
Since all the variables are in front of the matrices, they are treated as scalars. These are combined 
to form the familiar nondimensional parameters defined in Eqs. (3.21)-(3.24). Using the new 
matrices, vectors, and nondimensional parameters previously defined, Eq. (3.102) in 
nondimensional form becomes 

(M + fJ^°q + 2 p m nV V^q + (k + p m n 2 VV^q= - PfV (3.105) 


where 


in- ^ ext t 2 - F ext 

^ pL 3 " pvlL (3.106) 
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is a new nondimensional force parameter representing a generic nondimensional external force. 
Equation (3.105) next is transformed into modal coordinates, making it easier to place into state 
space domain. 

Modal Solution 

Equation (3.105) is now transformed from the discrete physical beam nodal coordinates, q , 
to the beam’s modal coordinates. At this point, the actual boundary conditions of the beam become 
integral. The equation, as it stands, is valid for any boundary condition; however, depending on 
the modes used, the equation is made specific for the boundary condition being examined. The 
exact equations for the simply- supported and the free-free systems are shown in Sections 3.2.5 
and 3.2.6, respectively. 

In Section 3.1, a modal transformation was explained for the continuous formulation. The 
form for the discrete formulation is the same. Instead of defining a continuous mode shape <pi(x), a 
discrete modal matrix 0, is defined. The actual modal coordinate rjj, is the same whether it is 

defined by the continuous mode shape and the beam's continuous displacement field u , or by the 
discrete modal matrix and the beam’s discrete displacement field q. 

q= L(prf (3.107) 

A description of some characteristics of the modal matrix follows. The actual 
transformation from physical coordinates to modal coordinates is completed. The resulting 
equation is nondimensional, is in modal coordinates, and is easily placed into state space domain. 

Modal Matrix 

The transformation from generalized coordinates to modal coordinates for the continuous 
formulation was shown in Eq. (3.11). A similar transformation, shown in Eq. (3.107), is valid 
for nondimensional vectors. The mode shapes used are in the form of a modal matrix and, like the 
mode shapes already used, this modal matrix does not have any dimensions. Unlike the scalar 
operation shown in Eq. (3.11), Eq. (3.107) represents a matrix equation. The modal matrix 
consists of the eigenvectors of the simplified system 
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(3.108) 


M q + K q = 0 

Equation (3.108) describes the discrete motion of the flexible beam without the moving 
mass. The number of modes present corresponds to the number of the system’s degrees of 

freedom. 

Transformation of Eq. (3.105) 

Using the substitution shown in Eq. (3.107), Eq. (3.103) is transformed into modal 
coordinates. A discrete version of Galerkin’s method, outlined in Section 3.1, is used. The 
resulting equation is premultiplied by the transpose of the the modal matrix to reduce the modal 
reduction error. This process is the matrix equivalent of using the continuous mode orthogonality 
to drive the error of the approximation to zero. For clarity, a definition of certain relationships 

follows. 

First, the modal matrix # is orthogonal with respect to the mass matrix, and is also mass 
normalized. This state leads to the following two definitions: 

<p T M<p=E (3.109) 

<p T K <p = A (3.110) 

E has previously been defined in the nomenclature. A is the nondimensional matrix of eigenvalues 
corresponding to the system shown in Eq. (3.108). The A matrix is different depending on the 
boundary condition being examined. The dimension of the modal matrix is equivalent to the 
system’s degree of freedom. For example, a free-free beam that is divided into ten equal beam 
elements has twenty degrees of freedom. It is not numerically efficient to retain all of these 
modes; therefore, before doing any numerical evaluation, the modal matrix and the corresponding 
eigenvalue matrix are reduced to retain a small number of nodes. In the actual numerical analysis, 
three flexible modes are retained. The reduced matrices are identified by a subscript r. 

Second, the following relationships are defined for the weighting vectors with the reduced 
modal matrix. 


Vjs V T <pr 


(3.111) 
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(3.112) 


T T 

V T)tt= V 44 ^ (3.113) 

Third, the modal damping matrix is developed. Since the modal decomposition of the 
discrete system is equivalent to the modes used for the continuous system, the modal damping 
matrix is the same. The damping due to the motion of the beam is assumed to be diagonal where 
the elements are defined by Eq. (3.32). The modal damping matrix is rewritten below for 
convenience. 


C q -2 C(XArt 12 (3.114) 

The matrix of Eq. (3.114) is defined in terms of the reduced matrix containing the eigenvalues 
defined in Eq. (3.110). 

Using the nondimensional modal coordinate T\ and the relationships defined in Eqs. 
(3.109)-(3.1 14), the nondimensional equation in modal coordinates becomes 

(E + (i m V v Vj)V + (2 £(XA r ) m + 2 HmnVr, V*) V 

+ (H + tt»« 2 H =-VfV T1 (3.115) 

Equation (3.1 15) is the discrete counterpart of Eq. (3.27). It is a nondimensional matrix 
equation rather than a continuous integral equation like that shown in Section 3.1. Though Eq. 
(3.114) is a matrix equation, it is not in the typical state space form easiest for numerical 
evaluation. 

3.2.4 State Space Representation 

Even with a matrix equation, the first step in forming a state space representation is 
choosing the state variables. Once the state vector is formed, Eq. (3.115) is transformed into the 
state space domain. 
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State Vector 


Two states variables are defined for the continuous system: the modal displacements and 
the velocities. Since the discrete analysis is already in matrix form, the state vector is written as the 
combination of two vectors 


x 



o 17 

V. 


(3.116) 


There is a difference between the number of elements used to model the beam n, and the number of 
modes retained to model the displacement of the beam N. The modal vectors have N elements, 
whereas the q vector contains n values. Usually only three flexible modes are retained, whereas 
there might be 100 beam elements used to capture the beam's physical displacement. 

Matrix Representation 


The mass, stiffness, damping, and force matrices that are used to describe the system's 
states are determined. Similar to the matrices obtained for the continuous formulation, the 
following matrices are a combination of a constant matrix and a time-varying matrix. The time- 
varying components are a function of the discretizations vector V. 


M = E + fi m Vjj Vq 
C = 2C(XA r ) m + 2n m n V n v\ % 


K = IA + H m n 2 V n vJ^ 


(3.117) 

(3.118) 

(3.119) 


It is interesting to see the similarities between the matrices shown here and the matrices 
shown in Eqs. (3.34)-(3.41) The matrices in Eqs. (3. 1 17)-(3.1 19) are the discrete counterparts of 
the previously shown matrices. The constant components of Eqs. (3. 1 17)-(3.1 19) are the standard 
mass, damping, and stiffness matrices obtained when developing a finite-element model of a 
flexible beam. The time-varying components actually model the dynamics associated with the 
moving mass. 


43 



The matrix representing the total discrete modal force applied to the beam must be formed. 
The actual force applied depends on the physical system being modeled. For example, if the beam 
is considered to be simply-supported, then a gravitational field is included as part o t e 
environment. However, when the free-free beam is examined, no gravitational field m included to 
correctly model the space envimnment. Sections 3.2.5 and 3.2.6 examine the discrete force matnx 

for each system in detail. 

The mass, damping, stiffness, and force matrices are used in Eq. (3.32) to solve for the 
modal displacements and the modal velocities. Using Eq. (3.107), the nodal displacements and 

velocities are obtained. 

This concludes the derivation of the general discrete formulation. Sections 3.2.5 and 3.2.6 
examine two specific systems. Section 3.2.5 analyzes the inertially fixed system, which is 
modeled using a simply-supported beam. Section 3.2.6 examines the inertially free system, 
which is modeled as a free-free beam. The results for each model are presented in Chapter 4, 
which compares them to the results obtained from the continuous formulation of Section 3.1. 

3.2.5 Applications to a Simply-Supported Beam 

Unlike the continuous formulation examined in Section 3.1.6, the discrete formulation for 
the simply-supported beam is no easier to formulate than the free-free beam. The simply- 
supported beam is examined and presented first so that its discrete methodology can be vail. dated 
against well-known results. The inertially fixed system is also used to determine which shape 
function. Unear or cubic, accurately models the beam with the smallest number of finite elements. 

Simply-Supported Beam using the Linear Shape Function 

For a simply-supported beam modeled with n finite elements there are 2n degrees of 
freedom; there are n - 1 translational degrees of freedom and n + 1 rotational degrees of freedom. 
As stated earlier, the linear shape function only uses the beam's translational degrees of freedom; 
therefore, each element's rotational degrees of freedom can be ignored. The reduced mass and 
stiffness matrices are derived from the mass and stiffness matrices shown in Eqs. (3.59) and 
(3.77). Since the beam is simply-supported, the first and last nodes are constrained to zero 
translation, leading to a mass and stiffness matrix with n - 1 degrees of freedom. 
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Because Ihe shape function is linear, its second spatial derivative is zero. However, as 
shown in Figure 4, there is a difference in slope between the two neighboring finite elements. This 
difference in slope leads to an important component of the inetlial force created due to the motton 
of the mass, which, for the linear shape function, is not present. To account for this mental force, 
which is proportional to the beam's curvature, an artificial impulse force is added to the equanon. 
In order to model this force, an impulse force is calculated as soon as the mass moves to the next 

element. 



Figure 4. Difference in slope of two neighboring finite elements. 


Using the value of the resulting force, an 
is applied over the entire element. 


fye ~ 


equivalent constant force defined in Eq. (3.120) 


fy_At_ 

At e ( 3 . 120 ) 


where At e is the time required for the mass to travel over one element. The actual impulse force, 
fyAT, is defined as 


f y At = m m v m ^Bq 


( 3 . 121 ) 


The vector B, defined in Eq. (3.122), is a central finite-difference operator that acts on the two V 
vectors, which determine the translations at each node. 

B = {vh-2 Vj + Vjj (3.122) 
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Using Eqs. (3.120), (3.121), and (3.122), the final equivalent force vector due to the beam's 
curvature is 



This force is added to the force term shown in Eq. (3.96). 

Comparing Eq. (3. 123) with the term that appears in Eq. (3.99) it seems that the vector B 

T 

is the linear equivalent of However, B is actually the second-order finite-difference 

approximation to the second spatial derivative of the V vector. 

Simply-Supported Beam using the Cubic Shape Function 

When the cubic shape function is used, the component of the inertial force due to the 
beam’s curvature results from the second derivative of the V vector. For the simply-supported 
beam depicted in Figure 2, the cubic shape function Vj, has 2n elements. These 2n elements 
correspond to the 2n degrees of freedom of the simply-supported case when the rotational inertias 
are included. 

Modal Matrix for the Simply-Supported Beam 

The modal matrix is evaluated by finding the eigenstructure of the simplified system 
depicted in Eq. (3.108). The dynamics are dictated by the mass and stiffness matrices. For a 
simply-supported beam, the matrices are reduced to eliminate the constrained degrees of freedom. 
The beam's boundary conditions specify that the translations at each end are zero. To address this 
condition, the first and last rows and columns of the matrices are eliminated. 

For the linear shape function, the mass and stiffness matrices contain only translations, 
which represents an n - 1 degree-of-freedom system. The reduced mass and stiffness matrices are 
variations of the matrices defined by Eqs. (3.59) and (3.77). 

The cubic shape function requires that both the translation and the rotations are present. 
Therefore, the reduced matrices come from the matrices shown in Eqs. (3.70) and (3.77). Note 
that if the system is cantilevered, it is not possible to statically condense out the rotational degrees 
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of freedom because they must be present in order to be eliminated to satisfy the boundary 
conditions. 

When modeling the simply-supported environment it is essential to include a gravitational 
field, the only external force applied to the beam. The gravitational force is proportional to the 
moving mass. The discrete nondimensional form of this force is 

VF* = - HgV (3.124) 

where Fg is the nondimensional gravitational parameter defined by Eq. (3.50). This external force 
matrix is the same matrix F, used in Eq. (3.32). 

3.2.6 Applications to a Free-Free Beam 

There are three differences between the formulation of the free-free beam depicted in Figure 
3 and the simply-supported beam evaluated in Section 3.2.5: 

( 1 ) The degree of freedom. 

(2) The modes used. 

(3) The external force. 

It has been shown that the cubic shape function is more efficient than the linear 

interpolation function when using the simply- supported beam. Therefore, in analyzing the free- 
free beam, only the cubic shape function will be used. Consequently, the mass and stiffness 
matrices are variations of Eqs. (3.70) and (3.77). The boundary conditions of a free-free beam 
state that the shear and moment at each end must be zero. These conditions do not constrain any 
degree of freedom. Therefore, the free-free beam and, correspondingly, the mass and stiffness 

matrices have 2n+2 degrees of freedom. 

These mass and stiffness matrices are used in forming the system’s mode shapes and 
eigenvalues. Before reduction, the modal matrix and the matrix of eigenvalues contains 2n+2 
elements. The first two modes correspond to rigid body rotation and translation and have zero 
frequency. The other 2n modes represent the beam’s flexible motion. 
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The free-free beam is used to model an inertially free system. In the space environment 
there is no gravitational field. If there is no external force applied to the beam their would be no 
response from the motion of the mass, unless an initial disturbance is used. The details of this 
vibration were outlined previously in Section 3.1.6. 

3.3 DISCRETE FORMULATION FOR THE FREE-FREE 
BEAM WITH MULTIPOINT OF CONTACT 

In the analysis thus far, the model has been a flexible beam with a rigid body attached at 
one contact point to the beam. In both the continuous and discrete formulations, a free-free beam 
and a simply- supported beam were examined. The simply-supported beam is used as a testing 
board for the discrete method outlined. The inertially free system is used to try and model the 
space environment of the Space Station-Mobile Transporter system; however, it still only models 
the mass as a wheel moving over the beam rather than as the more realistic train moving along a 
track. The next system analyzed, a free-free beam with multipoint of contact, is used to consider 
the train/track aspect of the SS-MT system. 

3.3.1 Mathematical Model 

The mathematical model used here, as before, is essentially the inertially free flexible beam 
with a mass moving at a constant speed along its length. For this analysis, however, the mass is 
attached at two points of contact (see Figure 5). 



Figure 5. Mass attached at two-points of contact 
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This final model represents the physical aspect of the mobile transporter as a train rather than 
simply a wheel as presented in Sections 3.1 and 3.2. The validity of the discrete methodology for 
this model is shown in Chapter 4, using the formulations developed in Sections 3.1 and 3.2. 

For the train/track model, only a discrete formulation is considered. When the two points 
of contact become infinitely close, the equations developed in this section should converge to the 
discrete equations for the one-point-of-contact case. This is proven to be true in Chapter 4. 
Therefore, since this method converges to a method that is already proven to be valid, there is no 
need to compare the results obtained from the following formulation with those of a continuous 
formulation. 

In Chapter 4, only an inertially free system is examined with two points of contact. In the 
following section a general derivation is developed. For the applications to the free-free beam see 
Section 3.2.6. 

The only difference between the discrete system examined in Section 3.2 and the system 
analyzed here is the manner in which the inertial force due to the moving mass is applied to the 
beam. Because there are now two points of contact, there are correspondingly two continuous 
forces due to the mass acting on the beam. Both forces must be made discrete and must be 
incorporated into the beam's equation of motion. 

As in Section 3.2, the analysis starts by discretizing the load applied to the beam. This new 
load is then incorporated into the beam’s equation of motion to obtain the equation of motion for 
the entire system. This equation is made nondimensional and placed into state space form for 
numerical computation. 

3.3.2 Equations of Motion 

The equation of motion for the train/track system depicted in Figure 5 is obtained with a 
series of steps. First, the discrete equation representing the displacement of the flexible beam is 
determined. Next, the force due to the moving mass is examined. Using compatibility, these 
equations are combined to form the system's equation of motion. Because this system is very 
similar to that depicted in Section 3.2, the following formulation is abridged to avoid repetition. 
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Beam Equation 


The discrete equation of motion for a flexible beam, shown earlier in Eq. (3.81), is 
rewritten here for convenience. 

Mq + Kq = F t (3.125) 

where, for an inertially free system, the mass and stiffness matrices are defined by Eqs. (3.70) and 
(3.77), respectively. F t is the total discrete force vector and is a wmbination of the external forces 
applied to the beam and the inertial effects of the moving mass. F t is the discrete form of Ft- As 
stated in Section 3.2.7, there are no external forces applied to the free-free beam. Therefore, 
without the loss of generality, the total force applied to the beam is only due to the inertial effects of 

the moving mass. 

Load Modeling 


The total load exerted on the beam due to the moving mass is the sum of the forces exerted 
by the two points of contact 


F t =fcl+fc2 


(3.126) 


The total load applied to the beam, due to the inertial effects of the moving mass, is the 
same as in the one-point-of-contact case. However, the load is now divided between the two 
points of contact. For simplicity, each point is assumed to carry half of the total load; therefore, 
the magnitude of the force at each contact point is half as great as the force exerted at the one 
contact point in Section 3.2. The force contribution from one of the contact points is 

Qml (3.127) 


where q m l is the discrete displacement of the first contact point. The load at the second point of 
contact is determined in a similar fashion. Substituting the actual values for the load contributions 
into Eq. (3.126), the load exerted on the beam due to the moving mass becomes 


m m ( d 2 q m l \ + % ( jPqmA \ 

{ 2 [dt* Us 2 \ at 2 ) abs ) 


(3.128) 
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Discretization Vectors 


The force shown in Eq. (3. 128) must be placed into a discrete form. In the derivation it is 
assumed that the distance between the two contact points is a multiple of the elemental length of the 
finite beam elements: 


A=jl e (3.129) 

where j is any integer between 0 and n. This assumption simply makes the bookkeeping a lot 
easier. 


The cubic shape function outlined in Section 3.2.2 is used to distribute the forces. The 
location of the first point of contact, xj , is determined. From this location, the position of the 

second point, x 2, is determined by knowing the speed of the moving mass and the spacing 

between the two points of contact. 

Because the two locations are dependent on each other, it is only necessary to follow one of 
the points; x \ has been chosen to locate the position of the mass. The force contribution from this 
point is discretized using the shape function shown in Eq. (3.92) and rewritten here: 

V 3 = (i -n 2 + 2 i 3 )v,Ai-2 4 2 + 

+ (^ 2 ■ 2 {*) Vj +2 + (-| 2 + V M ( 3 , i 30 ) 

The load contribution due to the second point of contact is discretized using a shape 
function of the same form but with different nodes that are determined by using the relationship 
between the two points of contact. If i represents the location of the first point xj, then i' 
represents the location of the second pointy. The relationship between i and V is 

i' = i-2j (3.131) 

Next, the weighting vectors for the second point of contact are formed. The new vectors are 
denoted by a prime and are formed using i' rather than i. 
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(3.132) 


Using these new V' vectors, an equivalent cubic shape vector for the second point of contact is 
determined. 

Vj = (j - 3 i 2 + 2 $ 3 ) Vi + (4 - 2 i 2 + 

+ (j | 2 - 2 4 3 )v'i+2 + (- r’ + v 'it3 (3.133) 

The equations in Section 3.2 were derived using a generic shape function V. The new 
shape function developed here, has the same basic format as V. Therefore, the discrete 
development that was shown in Section 3.2 is valid for the new vector Vj'. The equations that 
will be developed, however, are different because the force applied to the beam now has two 
components rather than one. 

In the following formulation, a generic V' is used to represent the new shaping function. 
Though a cubic shape function was developed (Eq. (3.133)) and is used in the numerical 
evaluation, the following methodology is valid for any order shape function. 


Some relationships are defined between the displacements of the contact points, the shaping 
functions, and the beams displacement field. 


Qmi = V T q (3.134) 

Qm2 - V' T Q (3.135) 


Using Eqs. (3.134) and (3.135), the absolute accelerations of the two points of contact are 
determined and substituted into Eq. (3.128). The resulting force is inserted into Eq. (3.100), 
resulting in the discrete force vector that is applied to the appropriate nodes of the flexible beam: 


F t = -^\vV T q + 2(^) + 

V' T q + 2(^) V'V'lq + (-ff V’V'^q] 
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where V'^ and V'^ represent the derivatives with respect to £ of the new shape function. For the 
duration of this derivation, V'£ and are shown in their symbolic form; their actual values are 
available in Appendix C. Now that the force due to the moving mass has been placed in a discrete 
form, it can be inserted into the discrete equation of the flexible beam. 


System Equation of Motion 

The discrete equation of motion for the system is obtained by combining Eqs. (3. 136) and 
(3.125). As in the previous cases, the terms involving the beam deflection q, are shown on the 
left-hand side of the equation even though they appear due to the moving mass. 

^f + ^L(w T + V' V' T )\q + 2^-^[w\ + V'V'l)q 


+ [K + f-^f(vV T 44+ V 


' V 'U)o = ° 


(3.137) 


Equation (3.137) represents a discrete equation of motion in physical coordinates. When the two 
points of contact are infinitely close to each other, Eq. (3.135) is equivalent to Eq. (3.105). This 
equation must be placed into a nondimensional form in terms of the beam s modal coordinates. 

3.3.3 Nondimensional Equations of Motion 

Noticing the similarities between Eqs. (3.105) and (3.137), it is trivial to place the latter 
equation into nondimensional modal form (see Section 3.2.3). Therefore, to avoid repetition the 
derivation is not repeated here. The nondimensional form in modal coordinates of Eq. (3.135) is 


(e+^(v„vJ+ v, v;))°n + 

2CM n + 2^n(v,vl (+ V v ^ 
+ v, ry)„=o 


(3.138) 


where the V'r\ vectors are found using the relationships outlined in Eqs. (3.11 1) - (3.1 12) for the 
new shape function. 
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Equation (3.136) is the multipoint of contact counterpart of Eq. (3.115); it is in matrix form 
but not in state space domain. Equation (3.135) is placed into the state space form in preparation 

for numerical computation. 


3.3.4 State Space Representation 

The only differences between Eq. (3.136) and (3.115) are the additional terms to the time- 
varying mass, damping, and stiffness matrices due to the second point of contact. The procedure 
for transforming Eq. (3.138) into state space domain is equivalent to that outlined in Section 3.2.4. 


State Vector 

The two states of the system are chosen as the modal displacements and velocities. The 
state vector is equivalent to the one presented in Eq. (3.116). Because a discrete analysis is already 
in matrix form, the state vector is rewritten in a slightly different manner than in the case of when 
the variables were scalars. This was first shown in Eq. (3.1 16) and is rewritten here: 

x = [ 7? V ] (3.139) 


Matrix Representation 

The matrices representing the system are a combination of a constant matrix and a matrix 
that varies with time. The constant matrices represent the dynamics of the flexible structure. Since 
the flexible structure modeled in this analysis is equivalent to that modeled in Section 3.2, the 
constant matrices are identical to those outlined previously. The matrices dependent on time model 
the dynamic interaction of the moving mass with the flexible structure. This interaction is the 
difference between the case in Section 3.2 and this case. The total matrices used to formulate this 
system in the form shown in Eq. (3.32) are 


M = E + ^(V,v; + V„V5) 

(3.138) 

C = 2 C(* a)" 2 *2^j-n[v , V'Jj + V'n V 'l() 

(3.139) 

K = AA*^Y^{v,v T riff *v n 

(3.140) 
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The above matrix definitions are valid for any boundary conditions. For the numerical 
evaluation of this system, only the free-free beam is examined. For the alterations needed to 
specifically examine a free-free beam, see Section 3.2.6. 
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CHAPTER 4 
RESULTS 

Chapter 4 presents the results obtained using the analysis outlined in Chapter 3. Section 

4.1 describes the layout of the results as well as the parameters used to create them. Section 4.2 
displays and discusses the results. 

4.1 RESULTS ORGANIZATION 

Section 4.1 describes the organization of Section 4.2, which displays the results obtained 
from the different formulations of Chapter 3. Three different formulations were developed in 

Chapter 3: 

(1) Continuous. 

(2) Discrete - one point of contact. 

(3) Discrete - multipoint of contact. 

In addition to these three different formulations, two specific systems were examined: 

(1) Simply- supported beam. 

(2) Free-free beam. 

To ease the complexity of the next section it is important to understand how the simulations 
are organized, what they are trying to accomplish, and what parameters are used to create them. 
To aid in this understanding. Section 4.1 is divided into two subsections. Section 4.1.1 discusses 
the goals of Chapter 4 and the way they are obtained. Section 4.1.2 examines the nondimensional 
parameters used to set up the simulations. For a description and listing of the computer codes used 
to perform the simulations, see Appendix D. 

4.1.1 Goals of Chapter 4 

Chapter 4 has three main goals: 

(1) Validate the discrete formulation for the one- and multipoint-of-contact cases with a 
simply- supported and a free-free beam. 

(2) Determine which shape function, cubic or linear, best models the system s 
displacement. 
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(3) Perform several informative parametric studies. These studies show the effects of 
the nondimensional parameters, developed in Chapter 3, on the system’s dynamics. 

Section 4.2 contains fifteen pages of plots. Each page displays either two, four, or eight 
plots, depending on the specific study being run. Figures 6 through 1 1 represent the simply- 
supported system. Figures 6, 7, 8, 10, and 1 1 represent time history of the nondimensional 
midspan deflection. Figure 9 presents a profile of the entire beam for different values of 
nondimensional time. The displacements shown in Figure 9 are also nondimensional. For the 
simply-supported beam all the displacements are made nondimensional by u$, the maximum static 

deflection of the beam. 

Figures 12 through 20 display results for the free-free system. For each study performed 
there are two sets of plots. The first plots, Figures 12, 15, 17, and 19, display the time history of 
the nondimesional deflection at the beam's left tip. The second plots. Figures 13, 16, 18, and 20, 
display the time history of the nondimensional deflection at the beam’s right tip. Figure 14 
presents a profile of the entire beam's nondimensional deflection. First, the organization for the 
simply-supported beam is explained; then, the free-free system is discussed. 

Simply-Supported Beam 

The first study performed for the simply-supported beam determines the shape function that 
best models the beam's displacement while using the least amount of finite beam-elements. To 
accomplish this, a comparison using a different number of beam elements was made between a 
continuous formulation and both the linear and cubic shape functions. The results presented in 
Section 4.2. 1 illustrate that the cubic shape function is better suited to model the displacement; 
therefore, in the following discrete simulations only the cubic shape function is used. 

The second study compares the discrete and continuous formulations for different values of 
the speed parameter a (see Eq. (4.5)). This study effectively validates the discrete formulation for 
the simply-supported beam. Consequently, the following simulations are performed for the 
discrete cubic formulation only. 

The third study presents snapshots that show the profile of the entire beam for different 
time frames. The first set of snapshots models the beam as the mass is moving along the beam. 
The second set examines the beam in free vibration, after the mass has left the beam. 
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Finally, a parametric study is performed for the simply-supported system, which examines 
the result of including the inertial effects of the mass versus simply modeling it as a moving force. 
First, different runs are completed for a specific mass ratio with different speed parameters. Then, 
the speed parameter is specified and the value of the mass ratio is varied. 

Section 4.2 contains the results outlined. They provide the information needed to reach the 
three goals set for the simply-supported system. The same basic tests outlined above, with the 
same parameters, were completed for the free-free system. This compatibility facilitates the 
comparison of the two separate systems. 

Free-Free Beam 

The discrete methodology for the free-free, one-point-ofO-contact system is validated. To 
accomplish this, a comparison is made between the discrete formulation presented in Section 3.2 
and the continuous formulation presented in Section 3.1, for various speed parameters. Based on 
the results of the simply-supported beam, only the cubic shape function is used. The speed 
parameters for the free-free beam (see Eq. (4.6)) are set to closely resemble those of the 
simulations performed for the simply-supported beam. Unlike the simply- supported beam, 
however, all the displacements are made nondimensional by the beam's length. 

Studies using the discrete formulation are also performed. First, snapshots of the beam are 
displayed for different time frames, as outlined for the simply-supported beam. The first set of 
curves models the beam with the mass travelling along its length; the second set of curves displays 
the beam in free vibration without the mass. 

The next set of plots explore the results of ignoring the inertial effects of the mass. If the 
inertial effects are ignored in the simply-supported beam, the mass still creates a gravitational force 
on the beam. When there is no gravity, no gravitational force is applied. Therefore, the two 
formulations compared in this study are: 

(1) With the moving mass. 

(2) Without the moving mass. 

The next study examines when it becomes important to include the inertial effects of the 
moving mass by varying the mass ratio for one speed parameter. Once again, the values of the 
parameters used are the same as those used for the simply-supported beam. 
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The final set of curves examines the multipoint-of-contact formulation. Different contact 
spacing simulations were compared, again using the same speed parameters that were used for the 
simply-supported beam. It is important to show that as the contact spacing approaches zero, the 
simulations approach the one-point-of-contact case. It is also interesting to see how the speed 
parameter alters the effects of the contact spacing. 

The above-mentioned curves help to achieve the goals that the simply-supported system 
could not obtain. Using the results from both the simply-supported and the free-free systems, the 
best shape function is determined, the discrete methodology is completely validated, and important 

parametric studies are performed. 

In the following simulations, the continuous formulation that was outlined in Section 3.1 is 
used for the validation of the discrete formulation. As an added assurance, an exact formulation 
developed by Kurihara and Shimogo (Ref. [5]) for the simply-supported beam was compared to 
the discrete simulation. The results obtained completely agreed with those obtained using the 

formulation detailed in Section 3.2. 

4.1.2 Parameter Discussion 

Each system is described by a stiffness parameter and the mass ratio between the structure 
and the moving mass. For both the simply-supported and the free-free systems, the stiffness of 
the system is represented by the nondimensional parameter A, and the mass ratio is represented by 
Mm . The external load for the simply-supported system is characterized by Mg- There is no 
external load applied to the free-free system. Before specifying the values of these parameters, it is 
important to define the speed parameter that was discussed in Section 4.1.1. 

Speed Parameter 

A convenient way to describe the flexible-structure/moving-mass system is to define a 
speed parameter a. 

T p T p Vpt 

a = 2t r = 2 L (4.1) 
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where Tp is the fundamental period of the system and •' ' £ was previously defined as the time 
required for the mass to travel the beam's length. To gain a more physical understanding of this 
parameter, a relationship between the fundamental period and the natural frequency of the beam is. 


which gives 


where 0)i is the beam’s fundamental natural frequency. 

Because a depends on both the speed of the moving mass and the frequency of the system, 
there are two different ways to look at the meaning of a. A low a represents either a system in 
which the mass travels at a very slow speed or a very stiff system. Conversely, a high a 
represents either a system where the mass travels at very high speeds or a flexible system. The 
actual physics of the system are the same regardless of how the speed parameter is interpreted. In 
this analysis, a is referred to as the speed parameter and consequently is used to describe the 

relative speed of the moving mass. 

The nondimensional natural frequencies of the simply-supported and the free-free beams 
can be expressed in terms of the nondimensional stiffness parameter. 


= k 2 U. 

(4.4) 

Q\ f = 22.4 fX 

(4.5) 


where Qj and ^.represent the natural frequency of the simply- supported and the free-free 
S J 

beams, respectively. 

Using Eqs. (4.3) and (4.4), a unique relationship is determined between the speed and 
stiffness parameters for the simply-supported and the free-free systems. The two speed 
parameters, oc s and <Xf , one for each specific system, are defined as 
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(4.6) 


a. 




«/= — ^ 7 = 

22.4 VA 


(4.7) 


Using Eqs. (4.6) and (4.7), the nondimen sional stiffness parameter A, is specified to achieve 
different speed parameters. Table 2 outlines the speed parameters used and the resulting value of 
A. for both beam's. Either the speed parameter a or the stiffness parameter A can be used to 
identify a specific case, because of the unique relationship between the two parameters. 


Table 2. Stiffness and Speed Parameters Used in Simulations. 


| Simply-Supported Beam 

Free-Free Beam 

a s 

A 

af 

A 

0.01 

1013.72 

0.1 

1.967 

0.2 

2.53 


0.492 

0.3 

1.12 


0.219 

0.4 

0.633 

0.4 

0.123 

0.6 

0.281 

0.6 

0.055 

0.8 

0.158 

0.8 

0.031 

1.0 

0.101 

1.0 

0.019 

1.2 

0.070 

1.2 

0.014 

1.4 

0.051 

1.4 

0.010 

1.6 

0.039 

1.6 

0.007 

2.0 

0.025 

2.0 

0.005 

3.0 

0.011 | 

3.0 

0.002 

4.0 

0.006 

4.0 

0.001 
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Remaining Parameters for the Simulations 

In addition to the speed parameter a (or stiffness parameter A), there are three other 

parameters of interest to the simulation: 

(1) The mass ratio, fi m - 

(2) The load parameter, fJ.g. 

(3) The contact point spacing, s. 

The mass ratio parameter p m , is defined the same for both the simply- supported and the 
free-free systems. The load parameter fig, is used only for the simply-supported system, and is 
not really an independent parameter, as will be seen. The contact point spacing parameter 5, is 
used only for the free-free system with multipoint of contact. 

As defined previously, p m represents the ratio of the moving mass to the flexible structure. 

This ratio is varied only in one set of simulations for each system (see Figures 11, 17, and 18). 
The purpose of varying p m is to determine the lowest mass ratio permissible to treat the mass as a 

moving force rather than as a moving mass. The value of the mass ratios vary from 0.01 to 2.0 
and are indicated in the legends of the appropriate curves. For the remaining curves, the mass ratio 
is kept at a value of 0.5; therefore, the moving mass is half as massive as the flexible structure. To 
model the mass as a moving force, the mass ratio is set to 0. 

In Section 3.1.5, a nondimen sional load parameter, fig, was developed to characterize the 
load applied to the system due to the gravitational force of the moving mass. The definition of Pg 
is rewritten from Eq. (3.51) as, 


_ rn m g 

p v2 (3.51) 


This load parameter in Eq. (3.30) would lead to nondimensional deflections u/L in the 
beam, for the simply-supported beam, it is more convenient to express the deflections in reference 
to the maximum static deflection of the beam u s . For the load m m g acting at the midspan, the 

maximum static deflection occurs at the midspan and is given as 
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Hence, to provide the nondimensional deflections u/us rather than the u/L would requre one 
to pick the new load parameter mg in Eqs. (3.30) and (3.50) as 





™ m g 
P v& 


m m g L 2 
48 El 


= 48 A 


(4.9) 


Using Eq. (4.9), the value of the load parameter pg is directly specified for the appropriate value of 
A or a specified in Table 2, or Eq. (4.6). 


The final parameter to be explained is s, the contact spacing parameter: 


100 


(4.10) 


where j is any integer. The parameter s represents the percentage of beam length by which the two 
points are separated. For example, if s is set equal to zero, only one point of contact is achieved. 
For s equal to 0.02, the two points are separated by a distance that is 2 percent of the beam length. 

Four different values of s were compared for different speed parameters: 


5 = 0.000 
s = 0.005 
s = 0.010 

5 = 0.020 


The parameters needed to describe the system have now been thoroughly explained. 
However, there are two more parameters that do not define the physical property of the system but 
do appear in Section 4.2. First, the simulations were run in terms of the nondimensional time 
parameter t and were run up to a value of r. = 2. Asa reminder, when T is equal to one, the mass 
has travelled over the entire beam length. Second, the parameter n defines the number of finite 
beam elements used. In Figure 6, n is varied. For the remaining discrete simulations, 40 beam 
elements are used. Finally, in the simulations three flexible modes were retained. 
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4.2 RESULTS DISCUSSION 


Section 4.2 presents the results for both the simply-supported and the free-free systems. 
Each graph is discussed as it appears in the text. 

The results that were outlined in Section 4.1 and their significance follow. Section 4.2 
contains two subsections, Section 4.2.1 discusses the simply-supported system, and Section 4.2.2 
expounds on the free-free system. The organization of each section was developed in the previous 

section. 

4.2.1 Results for the Simply-Supported Beam 

The plots are in terms of the nondimensional deflections ulus, and the nondimensional time 
parameter T = v m f/L. 

Figure 6 

Figure 6 compares the discrete formulation with both a cubic and a linear shape function 
and with the continuous formulation. Each curve was simulated for a equal to 1 .0. 

These curves are used to determine the shape function that accurately models the 
displacement with the least number of elements. Figure 6 (a) shows that for as little as 10 finite 
beam elements, the cubic shape function is almost identical to the continuous formulation. 
Therefore, for the rest of the discrete simulations the cubic shape function is used. It is also 
important to note that for sixty finite elements, almost no difference is seen in the three curves 

portrayed. 

Figure 7 

Figure 7 compares the discrete and continuous formulations for different speed parameters. 
These curves provide two useful pieces of information: they validate the discrete methodology for 
different speed parameters and they show the effect of varying the speed parameter. 

For each plot in Figure 7, the discrete and continuous formulations are almost identical. 
Figure 7 proves the validity of the discrete methodology for the simply- supported system. The 
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SIMPLY-SUPPORTED BEAM 






Figure 6. Discrete vs continuous for a = 1.0, /xm = 0.5. 
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Figure 7. Discrete vs continuous for varying a, fim = 0.5 
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plots were run up to a value of T = 2 in order to check the formulation when the mass is on the 
beam, and the free vibration when the mass leaves the beam. 

The speed parameters chosen to model the system range from a very slow-speed system of 
0.01 to a high-speed system of 1.4. Figure 7(a), a = 0.01, represents a system where the mass is 
travelling at a very slow speed. In this case, the beam sees the mass as a static force. The speed 
parameter is increased until a =1.4. By scanning the deflections as the speed parameter is 
increased, the effects of a on the system dynamics becomes obvious. 

For values of a less than one, the maximum effects of the moving mass occur while the 
mass is still on the beam. Even for slow speeds, i.e., a = 0.2, the motion of the mass affects the 
system dynamics, as indicated by the vibration that occurs even after the mass has left the beam. 
The largest deflection is detected when a = 0.6 (Figure (6(d)). For higher values of a, the 
damping terms due to the inertial effects decrease the midspan deflection. 

For speed parameters greater than one, the maximum effects of the moving mass occur 
after the mass has already left the beam. 


Figure 8 

Figure 8 displays the midspan deflection of the system for very large values of a. These a 
values represent the speeds that may be seen by a high-speed ground transport vehicle. As a is 
increased, the beam does not see the effect of the moving mass until values near t. = 1.5. It is 
suspected that as the speed parameter gets extremey large, the mass will have very little effect on 
the beam. This trend can be seen in Figures 8(a) through 8(d), especially in Figures 8(c) and 8(d). 


Figure 9 

Figure 9 displays two sets of curves. Each curve represents a profile of the entire beam at 
different time frames. For a speed parameter of 0.6, Figure 9(a) displays how the beam's 
deflection changes as the mass travels along its length. The maximum deflection occurs near 
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Figure 8. Very large values of a, /Am = 0. 
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Figure 9. Snapshot of beam with a = 0.6, fim = 0.5. 
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t- 1. Figure 9(b) simply shows the beam in free vibration after the mass has left the beam. As 
shown in Figure 8, for a = 0.6, no higher frequencies are present in the beam's vibration. 

Figure 10 

Figure 10 represents the first set of curves in a parametric study that examines the effects of 
the mass ratio parameter. In Figure 10, the mass ratio is set equal to 0.5, four different speed 
parameters are used, and two different formulations are displayed. The solid line represents the 
formulation where the inertial effects of the mass are included. The dotted line shows the 
deflection when the mass is treated as a moving force. 

For the static case, a = 0.01 (Figure 10(a)), no difference is detected between the two 
formulations. For the other speed parameters in Figures 10(b) through 10(d), however, a large 
difference in the two curves can be seen. The added inertial effects increase the effective force of 
the moving mass and consequently increase the maximum deflection of the midspan displacement. 

Figure 11 


Figure 10 examined the difference between a moving mass and a moving force formulation 
for a n m of 0.5. For all the speeds except the static case of a = 0.01 there is a difference in the 
two formulations. In Figure 11, however, the speed parameter is kept constant, a = 0.3, but the 
mass ratio fi m> is varied. 

This study determines when it is permissible to treat the travelling load as a moving force 
rather than as a moving mass. For a mass ratio as small as 0.01, there is a difference between the 
two curves. A significant difference, however, is not seen until the ratio reaches 0.05. The curves 
diverge as t approaches 1. 

In Figures 1 1(a) through 1 1(d), the mass ratio is still fairly insignificant (less than or equal 
to 0.1). Because the inertial effects add damping as well as an additional stiffness for small values 
the deflection shown by the moving mass simulation is smaller than the deflection predicted 
by the moving force simulation. But, as can be seen in Figure 10(e) to 10(h), for larger values of 
fifn the deflection predicted by the the moving mass formulation is larger than the deflection 
predicted by the moving force simulation. 
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SIMPLY-SUPPORTED BEAM 






Figure 10. MM vs MF for different a, (im = 0.5. 
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Figure 1 1. MM vs MF lor different /im. a=0.3. 
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For small mass ratios, the moving force assumption could be treated as a conservative 
prediction. However, for an accurate solution, the inertial effects should be included in the 
derivation. 

4.2.2 Free-Free Beam Results 

Each study that follows, except for Figure 14, contains two figures, one for each beam tip. 
The two corresponding figures are discussed simultaneously. The plots are in terms of the 
nondimensional displacement uJL, and the nondimensional time parameter t = v m r/L. 

Figures 12 and 13 

Figures 12 and 13 compare the discrete (cubic) formulation and the continuous formulation 
for different speed parameters. In every simulation the formulations are identical. These curves 
successfully validate the discrete methodology for the free-free system. 

Unlike Figure 6, the smallest speed parameter, a = 0.1, displays the largest amount of 
high-frequency vibration. To understand this phenomenon it is important to remember that the 
free-free beam was given an initial vibration prior to the release of the moving mass. Therefore, in 
reality. Figures 12(a) and 13(a) display that for the smaller speed parameters, the moving mass 
does not affect the system dynamics. As the speed parameter increases, the beam's deflection is 
increased but the beam's initial vibration is damped out. 

Figure 14 

Figure 14 is the free-free system equivalent of Figure 9. The beam's profile for different 
time frames is displayed. Once again. Figure 14(a) represents the beam as the mass travels along 
its length and shows how the beam is vibrating due to the initial kick it received prior to the 
presence of the moving mass. The top curve, t = 1.0, shows less vibration than the other 
intermediate curves. Figure 14(b) displays the free vibration of the beam after the mass has left. 
Since the beam is inertially free there is an absolute vertical motion of the entire beam. Figure 
14(b) also displays an oscillatory motion of the beam. 
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Figure 12. Disc, vs cont. for varying a, /im = 0.5, x = 0. 
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Figure 13. Disc, vs cont. for varying a, fim = 0.5, x = L 
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Figure 14. Snapshot of beam with a = 0.6, fim = 0.5. 
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Figures 15 and 16 


Figures 15 and 16 compare the simulations with the moving mass to the simulations 
without the moving mass. The comparison uses a mass ratio of fi/j j = 0.5 and different speed 
parameters. The simulations show a large discrepancy between the two formulations. Without the 
inertial effects, the beam is unaware that a mass is travelling across its length and acts like it is in 
free vibration. As a result, there is no rigid body translation or rotation. 

Figures 17 and 18 

Figures 17 and 18 determine when it is permissible to ignore the inertial effects that cause 
the discrepancy in Figures 15 and 16. The two systems, with moving mass and without moving 
mass, are compared with a constant speed parameter, a = 0.3, and varying mass ratios. 

Similar to Figure 10, a small difference is detected even for a mass ratio as small as 0.01. 
It is not until Pm = 0.05, however, that the difference between the two formulations becomes 
significant. It is interesting to note that the difference between the two formulations is more 
dominant at the beam’s left tip. The right tip deflections, predicted by the two deflections, are 
similar until a mass ratio is increased to 0. 1. 

Based on these simulations, it is not apparent that a moving mass assumption may be 
construed as a conservative approach for the free-free system. As concluded for Figure 10, to 
obtain accurate results, all inertial effects of the moving mass should be included. 

Figures 19 and 20 

Figures 18 and 19 focus on the formulation developed in Section 3.3, Discrete Formulation 
for Free-Free Beam with Multipoint of Contact. First, the curves are used to validate the 
multipoint-of-contact formulation. Second, the effect of the speed parameter on the multipoint-of- 
contact simulations is examined. 

For each speed parameter, simulations are compared for different values of s , the contact 
spacing parameter defined in Eq. (4.9). The spacing between the two points of contact decreases 
as s decreases. In both Figures 19 and 20, as the value of s decreased, the curves approach the 
one-point-of-contact simulation, s = 0. This trend is expected and consequently validates the 
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Figure 15. With and without moving mass for varying a , /zm=0.5 ( x=0 


79 







FREE-FREE BEAM 



NONDIMENSIONAL TIME 





Figure 16. With and without moving mass for vaiying a, /xm=0.5, x=L 
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Figure 1 9. Different spacing for varying a, fim = 0.5, x = 0. 
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Figure 20. Different spacing for varying a, /xm = 0.5, x = L. 
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multipoint-of-contact formulation derived in Section 3.3. As expected, for small-speed systems, 
the predicted deflections increase as the separation between the two points increases. 

Figures 19 and 20 display another trend: as the speed parameter is increased, the separation 
between the two points of contact has little effect on the system dynamics. Also, the right tip sees 
more variation than the left tip, which can be contributed to the extra time that a part of the mass is 
present on the beam due to the additional point of contact. The right tip sees this effect more 
significantly because the mass travels left to right across the beam. 
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CHAPTER 5 
CONCLUSIONS 


Chapter 5 discusses the important topics and conclusions presented in this thesis. The 
three different systems used to develop the methodology are examined, and conclusions reached 
for each system are discussed. Also presented in Chapter 5 is an overview of the different steps 
taken to develop the discrete formulation. Finally, Chapter 5 contains suggestions for future 
research. 

5.1 SYSTEM MODELS 

The primary goal of this research was to develop a methodology that may be used for 
modeling the inertially free Space Station-Mobile Transporter system using a discrete 
methodology. In order to validate the discrete formulation, a continuous formulation was 
developed. Three different systems were used: 

(1) Inertially fixed system (simply-supported beam). 

(2) Inertially free system (free-free beam) with one point of contact. 

(3) Inertially free system (free-free beam)with multipoint of contact. 

Each system was progressively more accurate in creating a model that resembles the SS- 
MT system. 

5.1.1 System (1): Inertially Fixed System (Simply-Supported Beam) 

The inertially fixed system was modeled as a simply-supported beam and developed for 
methodology verification purposes. The inertially fixed model used in this analysis is the same 
model that would be used for a heavy load travelling over a bridge. Therefore, many papers are 
available that provide a simulation of this simply-supported system. The results could 
consequently be verified by mere comparison with previously published results. 

The simplicity of the simply-supported modes enables easy formulation of the equations of 
motion. Not only are the equations easily formulated, but they can be fully expanded to display 
certain properties of the system (i.e., symmetry). Also, for a specific time, the mass, stiffness, 
damping, and force matrices can be determined by hand calculations, which provides a fast check 
for the computer code. 
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In addition to method verification, the effects of the nondimensional parameters 
representing the moving mass speed and the moving mass/flexible structure mass ratio were also 
examined. The dynamics of the simply-supported system can be easily interpreted by the results of 
the simulation; the situation is a bit more involved for the free-free system. Thus, the simply- 
supported system helped to develop a firm understanding of the moving mass/flexible beam system 
and, once understood, the formulations can be extended for more advanced problems. 

5.1.2 System (2): Inertially Free System (Free-Free Beam) with One Point 

of Contact 

The SS-MT system is designed to orbit around the earth. Therefore, to eventually model 
that physical system, a free-free system would be needed. In this formulation a free-free beam is 
used to represent an inertially free system. When the mass was attached to the beam at only one 
point of contact, the method was validated by comparing continuous and discrete formulations. 
Similar to the simply- supported system, the effects of the nondimensional parameters, which alter 
the physical properties of the system, were examined. 

The inertially free system with one point of contact successfully validates the discrete 
methodology that was developed. This system was also used to develop an understanding of how 
the nondimensional parameters affect the dynamics of the free-free system. 

5.1.3 System (3): Inertially Free System (Free-Free Beam) with Two 

Points of Contact 


The iteration to the inertially free one-point-of-contact system increased the complexity of 
the system but captured the train/track aspect of the SS-MT system. In the final formulation, the 
mass is attached at two points to the flexible structure. 

A continuous formulation was not used to validate the discrete approach for this system. 
Instead, the discrete one-point-of-contact case was employed. A comparison was made between 
the two-point-of-contact and the one-point-of-contact simulations for different contact point 
spacing. As the spacing between the two points approached zero, the two-point-of-contact 
simulation approached the one-point-of-contact simulation. Because the multipoint-of-contact case 
approaches a simulation that has already been proved, this provides credence to the discrete 
methodology for the multipoint-of-contact system. 
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5.2 CONCLUSIONS 


Using the three models described above, the discrete methodology was developed and used 
to simulate the dynamics of a mass moving over a flexible inertially free and inertially fixed 
system. The method was successfully validated for each system. An extensive parametric study 
was then performed and provides substantial insight in understanding this class of moving-mass 
problems. 

Six steps were taken to develop the discrete methodology: 

(1) Develop the discrete equation for the flexible structure. 

(2) Develop the discrete equation for the moving mass. 

(3) Using compatibility and finite-element shape functions,combine the two 
discrete equations into a discrete system equation of motion. 

(4) Place the system equation into a nondimensional form. 

(5) Perform a modal reduction on the entire system. 

(6) Place the nondimensional discrete system modal equation into a state space 
formulation for easy computational evaluation. 

In Chapter 3, each of these steps was discussed in great detail for the three stages of the 
model. The thrust of the research focuses on performing steps (2) and (3) with great accuracy and 
efficiency. The continuous formulations that were developed for method verification were derived 
in an identical manner. 

The analysis presented in this thesis develops a discrete methodology that will form the 
basis for formulating the SS-MT simulation. The simplified system examined here, an inertially 
free (or fixed) beam with a mass travelling along its length, is also solved using a continuous 
formulation for methodology validation. 

The method presented here is specifically designed with the SS-MT system in mind. This 
analysis formulated the dynamics of the Flexible Structure/Moving Mass system into a discrete 
form that is conducive to efficient computational analysis. The discrete formulation, which is in 
terms of nondimensional parameters that describe the necessary physical properties of the system, 
was placed in a form suitable for numerical integration. Modal reduction was used for 
computational feasibility. 
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5.3 SUGGESTED FUTURE RESEARCH 


Using the methodology presented here, a simulation of the Space Station-Mobile 
Transporter can be developed. Since the dynamics of the system may be cast in a discrete 
representation that would be suitable for modal reduction, the approach can be employed in the 
case of large dynamical systems: SS-MT-Space Shuttle. Because these systems will interact, the 
presented approach can be used to develop a simulation for dynamic interaction studies. 

For example, one possible scenario could be to examine the stability of the entire system 
when the mobile transporter is travelling along the space station while the shuttle is activating its 
attitude control system. The interaction between the Space Shuttle's dynamics, the attitude control 
system, and the SS-MT dynamics is just an example of the potential used of the approach 

developed in this thesis. 
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NOMENCLATURE 


NOMENCLATURE FOR CHAPTERS 1 THROUGH 5 


B 

C 

Co 

Cvar 

E 

El 

m e 

F 

Fext 


Discrete second derivative operator 
Damping matrix of the entire system 
Constant damping matrix 
Time-varying damping matrix 
Identity matrix 

Bending stiffness of the beam 

Elemental bending stiffness 

Force matrix of the entire system 

Sum of all continuous forces acting on beam 


F ext : 

Fjn 

F m 

Ft 

Ft 

fcl 

fc2 

fext 

frfext 

fye 

H 

i 

r 

K 

K 

K 

Ko 

K\ar 

k 

K t 

L 

le 


Discrete external force vector 
Continuous force acting on the moving mass 

Discrete force acting on the moving mass 
Total force acting on the beam 

Total discrete force acting on the beam 

Force due to first point of contact between moving mass and the beam 

Force due to the second point of contact between moving mass and beam 

Continuous external force applied to beam 

Modal external force applied to beam 

Equivalent impulse force 

Nondimensional integral where i = 1,2,3 and 4 

Location of first point of contact 

Location of second point of contact 

Stiffness matrix of the entire system 

Beam's stiffness matrix 

Nondimensional stiffness matrix of beam 

Constant stiffness matrix 

Time-varying stiffness matrix 

Global stiffness matrix with translational and rotational degrees of freedom 
Finite-element stiffness matrix with only translational degrees of freedom 

Length of beam 

Length of each beam element 
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M : 
M_ : 
M : 
M 0 : 
M t : 
Mvar • 
m : 


m e ■ 

mi 

mm 

N 

n 

Q_ 

Q 

Qm 

Qml 

Qm2 

s 

Ti 

Tp 

u(x,t) 


Ui : 
: 


V 
Vi 
V'i 

V T) 
Vl 
V3 

V'3 


Mass matrix for the entire system 
Mass matrix of the beam 
Nondimensional mass matrix of the beam 
Constant mass matrix 

Lumped mass matrix with only translational degrees of freedom 
Time-varying mass matrix 

Finite-element mass matrix with translational and rotational 

degrees of freedom 

Mass of each beam element 

Mass contribution of each node 

Mass of moving load 

Number of modes used to approximate the deformation 

Number of elements used to model the beam 

Beam's discrete displacement field 

Nondimensional discrete displacement field 

Discrete mass displacement field 

Discrete displacement field for first point of contact 

Discrete displacement field for second point of contact 

Contact point spacing in percent of beam length 

Trial function for cubic shape function where i = 1, 2, 5, and 4 

Fundamental period of the beam 

Continuous displacement field of beam with respect to an 

inertial reference frame 

Finite-element trial functions i = 1, 2, 3, and 4 

Displacement field of moving mass with respect to an inertial reference frame 

Mm = u m( x m(04) 

Maximium static deflection of the simply-supported beam 
Generic shape function 
Discretization vector i = 1,2,3, and 4 

Discretization vector for the second point of contact 
Modal form of V 
Linear shape function 
Cubic shape function 

Cubic shape function for the second point of contact 
Derivative of V with respect to E, 
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VvS • 

v m 

v m : 

x : 

xi : 
Xm(t) • 

Greek Letters 

a : 

«/ : 
a s : 

Pi : 
S( ) : 
A : 

: 

A : 
A : 
,V : 
V : 
m(t) : 
Vi(t) : 
<t> : 

4>r : 

<l>iM : 

(/,’ : 

t" ■ 
: 

P : 
01 : 


Second derivative of V with respect to ^ 

Modal form of 
Modal form of 

Relative velocity of the moving mass with respect to the beam 
State vector 

Position of the i th beam element 
Position along the beam of the moving load 


Speed parameter 

Speed parameter for the free-free beam 
Speed parameter for the simply-supported beam 
Parameter used for the mode shape of a free-free beam 
Dirac delta function 

Distance between the two points of contact 

Time required for the mass totravel over one finite beam element 

Nondimensional stiffness parameter 

Diagonal matrix of eigenvalues 

Reduced diagonal matrix of eigenvalues 

Vector of nondimensional modal displacements 

Modal displacement of the i^ 1 mode 

Nondimensional modal displacement of the i th mode 

Modal matrix 

Reduced modal matrix 

Mode shape of the M mode 

d<p 

dx 

d^± 

dx 2 

dU 

d x 4 

Mass per unit length of the beam 

Parameter used for the mode shape of a free free beam 
Nondimensional time parameter 
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Nondimensional gravitational load parameter 
Nondimensional force parameter 
Nondimensional discrete force parameter 
Nondimensional mass parameter 
Nondimensional frequency of the i*" mode 
Nondimensional frequency of the mode for the free-free beam 

Nondimensional frequency of the i th mode for the simply-supported beam 

Dimensional frequency of the mode 
Natural frequency of the beam 
Damping Ratio (C = -Ol) 

NOMENCLATURE FOR APPENDICES 

I c : Total mass moment of inertial of the system 

L : Lagrangian 

Mt : Mass of the beam and the moving mass system 

p(x,t): Vector locating deformed position of beam with respect to the 

undeformed position 

p m : Vector locating deformed position of moving mass with respect to the 

undeformed position 
Pm ~ Pm( x m( x ’ t )’ t ) 

r(t) : Vector representing the translation of the embedded body reference frame 

with respect to the inertial reference frame. 

Sc : Total static imbalance of the system 

T : Total kinetic energy 

Tb : Kinetic energy of the beam 

T m : Kinetic energy of the moving mass 

Vb : Potential energy of the beam 

v b : Absolute velocity of the beam 

v c : Absolute velocity of the moving mass 

v m : Relative velocity of the moving mass with respect to the embedded reference 

frame 

w(x,t): Vector representing the total displacement of the beam with respect to the 
embedded reference frame 


Pg 

Pfi 

Pf 

Pm 

Qif 

Qi s 

OH 

(01 

C 
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% : 


e 

a 


Vector representing the total displacement of the mass with respect 
to the embedded reference frame 
wm = w m( x m ( x ’ t )> t ) 

Angle representing rigid rotation between the embedded body frame 
and the inertial reference frame 

Skew angular velocity of the embedded frame with respect to the inertial 
reference frame 


co = 


0 -9 
0 0 
LO 0 


0 

0 

0. 
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APPENDIX A 

LAGRANGE FORMULATION 
OF CONTINUOUS FREE-FREE BEAM 

Appendix A formulates the continuous equation of motion for the free-free beam using a 
Lagrangian approach. The continuous formulation used to numerically simulate the motion of the 
free-free system presented in Section 3.1 developed the equations using a Newtonian approach. 

A Lagrangian formulation uses an energy approach whereas a Newtonian formulation uses 
force balance. If the individual forces of the system are known, it is easy to develop the equations 
using Newton's Laws of Motion. However, if the forces are difficult to identify, it is much harder 
to correctly develop the equations using a Newtonian approach. A Lagrangian formulation, on the 
other hand, uses the energy of the system. It is more difficult to formulate the equations but, if no 
algebraic errors are made, the equations are guaranteed to be correct. For this reason, for complex 
problems a Lagrangian formulation is developed to check the results obtained by the Newtonian 

method. 

Appendix A is organized in the same manner as the individual sections of Chapter 3. First, 
the equation of motion is derived. Next, the equation is made nondimensional and placed into a 
nondimensional form. This equation is then transformed from an equation in terms of the beam's 
natural coordinates, to an equation in terms of the beam's modal coordinates. Finally, a state space 
formulation is developed for this derivation. 

A.l MATHEMATICAL MODEL 

For this formulation a different model is used for the free-free system (See Figure A-l). A 
reference frame is embedded at the beam’s endpoints. With respect to this reference frame, the 
beam is simply supported at both ends. The reference frame is considered to undergo rigid body 
translation and rotation with respect to an inertial reference frame. 


In Chapter 3, the beam was able to move in any direction with respect to the inertial 
reference frame. There is no "body frame" embedded in the flexible structure. In that analysis the 
two rigid body motions fall naturally from the modal analysis and are not considered separately. 
That formulation is easier and more exact for numerical computations; however, the derivation 
developed here provides a good check for the equations because it is more thorough and less error 

prone. 
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Y 


Inertial Reference Frame 

Figure A- 1 . Model used for Lagrangian Formulation. 

A. 2 EQUATIONS OF MOTION 

The governing equations of motion are derived using basic energy principles. The kinetic 
and potential energy of the entire system is developed. These expressions are substituted into 
relations developed using Hamilton's principle. The resulting equations are known as Lagrange’s 
equation. 

To formulate the equations, the Lagrangian of the system, L, is used. The Lagrangian is 
the difference between the system's kinetic and potential energies. Calculus of variations is used to 
determine a function such that the integral of the Lagrangian takes on a minimum value (Ref. [19]). 
The resulting formulation is known as Hamilton's principle and is represented in symbolic form as 

S\ Ldt = 0 

h (A.l) 

Energy of the System 

In order to obtain the system's Lagrangian, both the kinetic and potential energy of the 
system must be determined. The system's kinetic energy is the sum of the beam's kinetic energy 
and the moving mass' kinetic energy: 
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(A. 2) 


T = T b + T m 


In determining the beam's kinetic energy, the flexible beam is considered to consist of an 
infinite number of beam elements. The moving mass is viewed as a ngid body. 

T b = \ f pVfcVfcdx 

2 Jo (A.3) 

T m = 2 mm v c . (A.4) 


In order to combine Eqs. (A.3) and (A.4) into a common kinetic energy expression, the 
velocity of the moving mass is rewritten using the special property of the dirac delta function 
shown in Eq. (3.14). The velocity of the moving mass is expressed as 

v c (x m> t)=\ v fc (x, t) S(x-x m )dx 

)o (A.5) 


Note, as discussed previously, this expression's validity is due to the fact that at x m the moving 
mass is firmly attached to the beam, ensuring equivalent velocities. Using the above relation, the 
total system’s kinetic energy is reformulated as 



v b + m m v T b \ b S(x - x m ) 



(A.6) 


where the dependence on the independent values are omitted for brevity. 

Because the beam's position vector is expressed in a reference frame that is moving with 
respect to the inertial reference frame, the inertial derivative contains two variables. The first 
variable reflects how the beam's position changes in time with respect to the embedded reference 
frame. The second value determines how the motion of this embedded frame changes with time 
with respect to the inertial frame. The expression for vfj is 
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(A.7) 


3 r .dp 

Vb= a? + Ji + 


cap 


where the two vectors r and p are defined as 

r. Vector locating the origin of the embedded reference frame 
with respect to the inertial reference frame. r = [r x ry r z ] T 
p: Vector locating the beam's deformed position with 

respect to the undeformed position. It is expressed with 
respect to the embedded reference frame. p = [xw 0] T 


There are two other vectors necessary to form the kinetic and potential energy expressions. 
They are: 


O'. Angle representing the rigid rotation of the embedded reference 
frame about the inertial reference frame. 
p m : Vector locating the moving mass' deformed position with 

respect to the undeformed position. It is expressed with 
respect to the embedded reference frame. p m = [x m w m 0] T 

Using these vectors and Eq. (A.7) the velocity expressions in the x and y directions for the 
beam are: 


(vik = r x - (A. 8) 

(\ b )y = r y +W+ Ox (A9) 

Since the reference frame embedded in the moving mass is moving, the velocity of the 
moving mass takes on a different form. It is expressed as 

(v c )x = ^ - 0w m + v m (A. 10) 
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(A. 11) 


( v c)y - r y 


+ w m + V 


dw„ 


dx„ 


+ 9 x„ 


where v m is the speed of the moving mass relative to the beam. Eqs. (A. 10) and (A.l 1) represent 
the moving mass velocity that would be obtained using Eq. (A.5). 


Next, the system's potential energy is defined. The system being examined is in a gravity- 
free environment. Thus, the only potential energy of the system is the strain energy due to the 
beam's deformation: 


rL 


V*= 


El 


( d 2 w 1 


dx 2 


dx 


(A. 12) 


where El is the effective bending stiffness of the beam. Note that the beam’s material properties 
are assumed to be homogeneous. The potential and kinetic energy expressions are used to form 
the Lagrangian, which is used in Hamilton’s Principle. 


Hamilton's Principle 


Once the system's potential and kinetic energy are known it is trivial (but tedious) to apply 
Hamilton’s principle. First, the Lagrangian of the system is formed. Substituting Eqs. (A.6) and 
(A. 12) into the Lagrangian expression yields Eq. (A. 13): 




L = 


Jo 


I 


1 

2 


v£ Vb + 


v[ v fc S(x - x m ) 



dx 


(A. 13) 


Equation (A.13) is substituted into Eq. (A.l). The variation of L is taken with respect to 
each of its dependent variables. In this system, L is a function of eight variables: 

L = 4. r y> ' r v ™m, Xm ) (A. 14) 
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where 0 and ( ), represent derivatives with respect to time and spatial position, respectively. The 
r and w vectors have been previously defined. The other variables are: 

It is necessary to take the first variation of L with respect to each of the variables shown in 
Eq. (A. 14). In symbolic terms Eq. (A.l) becomes 



where the actual expressions for the kinetic and potential energy, Eqs. (A.6) and (A. 12), have been 
omitted for clarity. It is desirable to have Eq. (A. 1 5) in a form where the only variations are of the 
actual variables, not their respective derivatives. This form is: 

f [(• • •) Sr x + (•••) Sr y + (•••) Sw + (■ • •) $6>] dt -0 

it, (A. 16) 


To obtain the form outlined in Eq. (A. 16), each term in Eq. (A. 14) is integrated by parts. 
This type of integration separates two functions. One function becomes a differential and the other 
is integrated upon. The form of this type of integration is: 


1 


udv = uv - 


I 


v du 


(A. 17) 


The function chosen to be u takes on a differential form in the new integral. The function chosen 
to be dv is in its integrated form, v, for the new integral. The term in front of the new integral is 
evaluated at the endpoints of the integral. 


This type of integration is performed on the terms in Eq. (A. 14). The variation is chosen to 
be the dv function and the corresponding differential is the u function. The terms that are evaluated 
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at the endpoints of the integral are the beam’s boundary conditions. After the integration by parts, 
the equation is in the "strong form" (Ref. [18]). 



These four integrals are separated to obtain four equations of motion. 

Final Equations of Motion 

For Eq. (A. 19) to be true, each integral must vanish independently. To have each integral 
vanish for any arbitrary time period, the actual integrands of each integral must respectively go to 
zero. This leads to four equations. 

The four equations are obtained by following three steps. First, the actual kinetic and 
potential energies outlined in Eqs. (A.6), (A.7), (A.8), (A.9), (A.ll), and (A.. 12) are used in Eq. 
(A. 14). The resulting expression is integrated (by parts) to obtain the form shown in Eq. (A. 15). 
The additive property of this integration is used to separate the result into the four equations of 
motion, shown in Eqs. (A.20) - (A.23): 

(m m + pL) r x = 0 (A.20) 
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{m m +pL)r y +( 
2 dw 

W + m m \ 2 m + 


pL 2 


+ rn m x m 


2 m m v m 


<^r ohc 


6+2 nt m v m 



0 + 

w dx - 0 


(A.21) 


pL' 


+ /W/n 


■ {PL 3 

r y + rr 


dw 


+ m m xjh\ 6+2 m m \ m x m 6 + 

•L 


Xm VV + trim %m ''m ^ 2 TTim XmS , 


d?W 
dt dx 


} @Xm 


I 


+ I pxwdx -0 


(A. 22) 


pL r y + (pL 2 + m m x m )6 + 2 m m v m 6 + 


m m w + tn m vfn 


d 2 w 

'dx 2 


+ 2 m m v„ 


d 2 W ' 
dt dx ^ 


+ p w + El = 0 
'&>X m d X 


(A.23) 


Equations (A.20) and (A.21) represent Newton's second law in the horizontal and vertical 
direction, respectively. Equation (A.22) states that the sum of the moments around the origin of 
the embedded reference frame is zero. Equation (A.23) is the partial differential equation 
describing w(t), the lateral vibration of the beam. These equations could have been written 
directly using Newton's Law of motion (see Section 3.1); however, it is important to account for 
all of the forces. The energy approach might be more time-consuming than if Newton's Law were 
immediately applied, but if done carefully it assures that all forces acting on the system have been 
represented correctly. 

Equations (A.21) through (A.23) are next placed into the beam's modal coordinates. In 
this form it will be easier to compare the terms obtained in this derivation and the ones obtained 
during Section 3.1. Section 3.1 was derived for any boundary condition. This analysis is specific 
to the free-free beam. 


Modal Solution 

The beam's total deflection must be determined. The total deflection is the vector sum of 
the rigid body motions and the flexible motions of the beam. The horizontal motion of the system. 
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Eq. (A.20), is decoupled from the other three motions, thereby not playing a role in the total 
deflection. The other three equations are completely coupled. 

A linear superposition of modes is used to solve the three coupled equations, Eqs. (A.21), 
(A. 22), and (A.23). As in Chapter 3, Galerkin’s method is used to reduce the error of the 
approximation. However, unlike the lone partial differential equation shown in Eq. (3.10), there 
are two ordinary differential equations and one partial differential equation. A modal substitution 
for the beam’s vibration is used in all three equations, but Galerkin’s method is only applied to the 
partial differential equation. 

The substitution shown in Eq. (3.11) is used for the beam’s lateral vibration, w(x,t). The 
modes chosen must reflect the beam's lateral vibration only and not the total displacement as was 
modeled in Chapter 3. After the modes are chosen, the three equations are transformed into modal 
coordinates and Galerkin’s method is used where applicable. 


Modes Used 

The modes used are those of a simply- supported beam. This may seem incorrect since the 
beam itself is considered to be inertially free. However, this derivation was formulated so that the 
beam is simply supported with respect to the embedded reference frame and the embedded frame 
undergoes the rigid body motion with respect to the inertial frame. Therefore, the choice of 
simply-supported modes for the vibration of the beam is justified. The modes for the simply- 
supported beam are shown again for convenience. 

<t>i{x) = sin (a. 24 ) 

These modes are nondimensional and orthogonal; however, it is noted once again that they 
are not orthonormal. Using these modes, Eqs. (A.21), (A.22), and (A.23) become 

(m m + pL) r - y + J^- + m m x^0+ 2 m m v m 8 + p E rjij <pi dx + 

N \ 

+ m m E U rji + 2 v m <t>i T}i + v 2 m <pi t]i)@ *„= 0 (A.25) 
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+ ftl m X m Vy + 


plJ_ 

3 


+ m n xl 0 + 2 rn m v m x m 6 


N 

+ p E 

i = 1 



x fa dx 



(A.26) 
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f N N \ 

El E <pr rji + pr y + pxd+ p E <l>i jji 

i = 1 i = 1 > 


dx 


Jo 


+ m„ 


+ m m (t)j(r : 


+ x m d+2v m 6 


y T 


») 


(A. 27) 


<t>j 


N 

I E ( (pi Vi + 2 V m (pi Vi + v 2 m (pi Vi ) 


i = l 


k® X m 


0 


Equations (A.25), (A.26), and (A.27) represent the motion of the inertially free flexible 
beam with a mass moving along its length. The equations are in terms of the modal coordinates 
but are still in dimensional form. Therefore, the next step is to make them nondimensional. 

Some nondimensional integrals are needed to place the equations into a nondimensional 
form. Two of the nondimensional integrals were already defined in Eqs. (3.15) and (3.16). They 
are rewritten below for convenience. 


h{iJ) = f I <Pi <Pj dx 

Jo (A.28) 

h(iJ) = E 3 [ <pr (pjdx 

Jo (A. 29) 
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The two other integrals needed are actually just special cases of Eqs. (A.28) and (A.29) that 
incorporate the rigid body modes. Since these motions are not accounted for in the mode shapes, 
as in Chapter 3, the integrals must be defined separately. They are 

uP 

. 

II 

t' , b» 

^ £ 

&• 

(A. 30) 

U‘\ = -(■ f xipidx 

L 2 )o 

(A.31) 

Equation (A. 30) represents the rigid body translation and Eq. (A.31) is the rigid body rotation 
(i.e., <pj = 1 and <pj = x, respectively). The values of the /; and / 2 , for the simply-supported mode 
shapes, were found in Eqs. (3.45) and (3.47). These values together with I 3 and I 4 . given below 

ia>a=2 

(A.32) 


(A. 33) 

l (/) = 1 - cos±n_ 

W ITt 

(A. 34) 

U (i) = - £n ^ JjL 

w in 

(A. 35) 


For convenience the following variables are defined 


M t = pL + m m 

(A. 36) 

_ pL 2 

S c = — r — + mm Xm 

(A. 37) 

pL 3 , 

Ic “ j + 

(A. 38) 


1 1 1 



where M t is the mass of the beam and the moving mass, S c is the total static imbalance, and I c is 
the entire mass moment of inertia. Using the above identities and the four nondimensional 
integrals, the three equations of motion are rewritten in a further condensed form: 


N 

M,r y + S c d + 2m m v m d + pL E h(0 Vi 

i = 1 
N 

+ m m E ($ Vi + 2 v m (pi rji + w 2 m (pi Vil@ x m = 0 
i = 1 


(A.39) 


N 


S c r y + I c 6 + 2m m x m \ m 6 + p L 2 E h(0 Vi 


i = l 


+ Ttl m Xn 


E ((pi Vi + 2 v m (pi Vi + v 2 m (pi Vi\@ x m = 0 


i = l 


(A.40) 


(pL h{j) + rn m <p j@ J r y + ( pL 2 1 4 {i) + m m x m (p j@ J 6 + 

N N 

2 rn m v m <pj 6 + p L E I\(iJ) Vj + EL .X hUJ) Vj 


N 


+ m m | (pj E ((pi Vi + 2v m (pi Vi + vli to Vil 


© Xm 


= 0 


(A.41) 


It is important to note that the equations, as they stand, are valid for any mode shape used. 
The mode shapes should satisfy the geometric and force boundary conditions. 


This analysis uses modes shapes for a simply-supported beam, i.e., sin modes. With 
respect to the embedded frame, the beam is simply-supported; therefore, the modes do satisfy the 
geometric boundary conditions. They do not, however, satisfy the force boundary conditions for 
the actual free-floating beam. Using the simply-supported modes, there is a shear force at the ends 
of the beam. For a free-floating beam there are no moments or shears at the endpoints. The 
natural frequencies of the above modal system, with the rigid motions constrained, resemble the 
natural frequencies of a free-free beam, but the shapes of the sine modes do not accurately model 
the motion of the actual beam. It is better to use the free-free mode shapes of the beam as in 
Chapter 3, which satisfy all boundary conditions at the ends. 
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For these reasons, the above derivation is not used to numerically determine the beam's 
total deflection due to the motion of the moving mass. Even though this derivation is not used to 
compute the actual results, it is presented as a check to the derivation of Chapter 3. 

A. 3 NONDIMENSIONAL EQUATIONS OF MOTION 

Equations (A. 39), (A. 40), and (A.41) are placed into a standard nondimensional form 
using the procedure outlined in Chapter 3. The same nondimensional parameters that were defined 
in Eqs. (3.21)-(3.24) appear. This derivation is specific to the inertially free system; therefore, 
there is no external force applied to the beam. Since the mode shapes used are sine waves, the 
respective derivatives are easy to evaluate and have also been included. The final three 
nondimensional equations are: 


N 


(1 + Hm) °r° y + Ur + fJm f) 8 + 2 m m 8+ X I 3 (i) °ifi 

& i - 1 

+ m m X sin inx rji + 2 in cos inn rji - (inf sin inx 77 , j = 0 


(A.42) 


11 \ // \°° o 

(j + UmtjVy +(|- + /I m r 2 ) 8 + 2 fi m r 8+ X I 4 {i) °V°i 

[ N 

+ Li m t\ X sin inx 

Vi = 1 



(A.43) 


(h(j) + sinjnx) °r y + (I 4 (j) + m m xsin jnx) 8 

O N oo N 

+ 2 Hm Sinjnx 8 + X I\ (ij) °i?j + X X h(ij) Vj+ 

( N 

I . ° ° O 

/i m sin jnx X Sln lKX r l‘ + 2 in cos inx Vi - W sin inx rj; 

\i= 1 


for j=l,..,N 


= 0 


(A.44) 


These three equations are the nondimensional counterparts of Eqs. (A. 39), (A.40), and 
(A.41) and correspond to the expanded version of Eq. (3.27) Equations (A.42) and (A.43) 
correspond to the rigid body translation and rotation, respectively, and Equation (A.44) 
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corresponds to the flexible modes. Equation (A.44) along W ith r y = 0,6=0 comesponds ,o .he 
simply-supported beam case developed in Chapter 3. 


Though these equations will not be used for numerical analysis, they are still placed into a 
state space" form because i, is easier to see the similarities and differences between the denvatton 
presented here and the one reviewed in Chapter 3. 


A. 4 STATE SPACE REPRESENTATION 


The three equations above, Eqs. (A.42), (A.43), and (A.44), are a coupled set of N first- 
order equations. The three equations are combined into one matrix equation tons depends, 
the system states. In a state-determined system, all informatton needed about the system 
summarized in a finite set of variables. First, tire state variables are placed tnto a state vector. 
Using this vector, a matrix representation of the system is obtained. 


State Vector 

The following two vectors are defined. 

x l - { r y 6 t\j Vn ) (A.45) 

\T 

x2 = \fy 0 m-'-m (A. 46 ) 

The state vector x, is a combination of xl and x2 

* = IS) < A 47 > 


As stated previously, this state vector contains all the relevant information about the system. 

Matrix Representation 

Using the state vector of Eq. (A.47), the mass, damping, and stiffness matrices are 
derived. As indicated earlier, there is no external force applied to the beam; therefore, there is no 
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external force matrix. The mass matrix is the only symmetric matrix. All the matrices are time 
varying. 

For easy comparison to the earlier analysis, each matrix is divided into four different 
matrices. Three of the matrices have contributions from the rigid body modes. The fourth matrix 
contains contributions from the flexible modes only. These last matrices are identical to the 
matrices shown in Eqs. (3.47), (3.48), and (3.49). 


The mass matrix is: 


M = 


Ml 
. M2 t 


M2 
M4 . 


2+N x2+N 


(A.48) 


where the three matrices are defined as 


Ml = 


1 + Mm 


^-+Mm^ 


2x2 


(A.49) 


M2 = 


h(l) 

+ fi m sin nr 

U{1) 

+ (i m r sin nr 


h(N) 

+ fi m sin Nnr 

im 

+ fx m Nrsin Nnr 


(A. 50) 


M4 = 


fJ-m 


h (U) + 

sinnr sinnr 


sym 


Ii(N,l) + 

sinnr sinNnr 


NxN 


lj(N ,N) + 

fi m sinNnr sinNnr _ 


(A. 51) 
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Similarly, the damping and stiffness matrices are 


where 


c = \ Cl C2 ' 
C3 C4 . 


and 



K2 
K4 . 


(A.52) 


Cl = 


0 

0 


2 Mm 
2 Mm* . 


2*2 


(A.53) 


C2 = 


2 sin nx 
2 /i m r n sin nx 


2 Mm Nn sin Nnx 
2 Mm t Nn sin Nnx 


N x 2 


(A.54) 


C3 = 


0 

0 


2 Mm sin nx 


2 x N 


2 Mm sin A/ttt 


(A. 55) 


2^. m nsinnxcosnx 

+ 2 £ x 

2 


C4 = 


2fi m NnsinnxcosNnx 


\N x N 


2fi m nsinNnx cosnx 


2 fi m N nsinN xcosN nx 
+ 2 £ x 

t/XIi(NJ*)I 2 {NJ4) 

2 


(A. 56) 


where modal damping has been added to the fourth damping matrix. 
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The stiffness matrices are 


Kl = 


0 O' 
.0 0 . 


2x2 


(A. 57) 


K2 = 


- Hm (nf sinnx 

- fi m x ( jtfsinrcT 


- fj. m (NnfsinNnx 
- ( NnfsinNnx 


N x 2 


(A. 58) 


K3 = 


0 0 

o o 


2xN 


(A. 59) 


Ul{N,\y 

fi m (N nfsinKt sinNnx 


X Il{N,N)- 

fj. m (N K^sinN kx sinNnx _ N x 2 

(A.60) 

These matrices form the mass, stiffness, and damping matrices that may be used to 
simulate the motion of the moving mass system. The values of//, /2, and 14 for the simply- 
supported case are given by Eqs. (A.32)-(A.35). In order to obtain the total deflection, the ngid 
body motions must be added to the beam's flexible deformation. This analysis functioned well as 
a check for the Newtonian method and for that reason has been included here. This formulation 
was not used for any numerical evaluation. 


K4 = 


A/ 2 (U)- 

fj. m n 2 sinnx sinnx 


XI 2 { i,n)- 

fi m n 2 sinN nx sinnx 
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APPENDIX B 

RUNGE KUTTA INTEGRATION SCHEME 


The most important and significant problems in engineering are formulated in mathematical 
terms as a function that satisfies an equation containing the derivatives of the unknown function 
(Ref. [20]). Such an equation is termed an ordinary differential equation. The theory of 
differential equations dates back to the seventeenth century with the beginnings of calculus (Ref. 
[21]). Solutions for different types of differential equations were derived by such great 
mathematicians as Newton, the Bernoulli brothers, and Euler. The French mathematicians, 
Lagrange and Laplace, also made great contributions toward the solutions of ordinary differential 
equations (Ref. [20]). 


However, there are still a number of ordinary differential equations for which no analytical 
solution has been found. Instead, numerical integration is used to solve for the function satisfying 
the ordinary differential equation. 


Before computers were invented, numerical integration was hand-calculated. The most 
well known numerical integration scheme is Euler's method, which can be derived by forming a 
Taylor Series expansion of the ordinary differential equation. For example, consider: 


x = F[x,t) 


(B.l) 


This ordinary differential equation is the type that must be solved to simulate the motion of the 
flexible-beam/moving-mass system. 


A Taylor series expansion of Eq. (B.l) leads to 


Jjt+i =x k + F k Ar + H.O.T 

where the subscript indicates the time level: 

Xk'. x at time t 

*k+l‘ *attimer + A/ 

Fk- Function evaluated at xk and t 


(B.2) 
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Equation 


ion (B.2) is rearranged to obtain the forward Euler method of numerical integration. 

(B.3) 


F*oitol -St 


At 


Because only first-order terms are retained, this method is termed first-order accurate. Tbe tide 
accurate" implies that the accumcy is proposal to the * order of the ume step, A*. 

Using hand calculations, Eq. (B.3) can be iterated. Because i, is 
however, many iterations are needed to converge the method to dte correct 

throughout the ages, more advanced numerical integration schemes have been de 

new methods are more complicated and tend to have a higher order of accuracy. 

It would not be practical to employ these more advanced methods by hand; but, using a 
computer!" can be implemented very easily. As computing power increased, more an mom 
numerical integration schemes were developed. Three of these methods^ are dtscussed; the step- y- 
step, the predictor cotrector, and the alternating direction implicit methods. 

The step-by-step method examines the function at intermediate time steps. These 
intermediate values are weighted and combined to obtain die appropriate answer. The accuracy of 
die step-by-step method increases as the number of intermediate steps used increases Method 

Le known as explicit schemes fccause they use die value a, a previous time frame to 

determine the value at the current time. 

Conversely, an implicit scheme uses the value at the cinrent time to predict the value at a 
future time The combination of an implicit and an explicit scheme is used to form a predicmr- 
corrector method. An implicit scheme is used to predict die future value and then an explicit 

scheme is used to correct this predicted value (Ref. [22]). 

The alternating-direction implicit schemes are variations on the finite-element method. 
Unlike the first two methods, this scheme works best on a mesh of elements rather than on a string 
of nodal points. This type of numerical integration is only practical for two-dimensmnal (or^ 
problems. For half of the time step, a sweep is made in one direction, i.e„ r, or 
of the time step, a sweep is made in the other direction, i.e., y (Ref. [22]). 
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As well as the three categories discussed here, there are dozens of numerical integration 
schemes These three were chosen to provide an overview of the different integration methods 
available. However, to simulate the motion of the nexible-beam/moving-mass system, a step-by- 
step integration scheme, a Runge Kutta integration scheme, was chosen. 

Four intermediate steps were used so the scheme would be fourth-order accurate. The 
fourth-otder Runge Kutta integration scheme for the type of ordinary differential equation shown 

in Eq. (B.l)is (Ref. [19]): 

x k+i=x k + L(b l + lb 2 + 2b 3 + b 4 ) (B .4) 


where 

b x = AtFix k , t k ) (B.5) 

b 2 = At Fic (x* + £ b\ , t k + (B .6) 

b-i = At F k (x* + ^b 2 , t k + ( B .7) 

b 4 = AtF k (x k + b 2 ,t k + At) (B.8) 

This integration scheme was also used to numerically integrate the nondimensional integrals 
for the free-free mode shapes. The values were then used to formulate the constant mass, 
damping, and stiffness matrices, the values of which are shown in Appendix C. 

For numerical stability, the time step used in the integrations was smaller than one tenth of 
the shortest period for each system. If there are many modes retained, the shortest period can 
become very small. For computational efficiency only the first three flexible modes were retained 

for each system. 
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APPENDIX C 

NUMERICAL VALUES FOR THE FREE-FREE BEAM 


For clarity, the spatial derivatives of the shape functions, linear and cubic, used in the 
discrete formulation were kept in their symbolic form. The constant mass, damping, and stiffness 
matrices were also left in an unexpanded form because the nondimensional integrals needed to 
form these matrices were only obtained as numerical values as a result of a numerical integration 
scheme. This appendix presents the numerical values of these matrices and develops the 
appropriate derivatives of the shape functions with respect to the weighting function £,. 

C.l NUMERICAL VALUES OF THE CONSTANT MATRICES 

The mode shapes used to model the free-free beam consisted of trigonometric and 
hyperbolic trigonometric functions. Due to the complexity of these mode shapes, the values of the 
nondimensional integrals needed to form the constant mass, damping, and stiffness matrices were 
determined by using the Runge Kutta integration scheme outlined in Appendix B. 

As a refresher, the mode shapes used to model the free-free beam are rewritten: 

01 = 1 (C.l) 

02 = T- 1/2 (C.2) 

(pi = cos pi x + cosh pi x - Oi (sin Pi x + sink pi x) 3 < i < N (C.3) 

where, as before 

7 = i 

L (C.4) 

The numerical values of the parameters pi and oi are shown in Table C-l for the first three 
flexible modes (Ref. [17]). 


123 


PRECEDING PAGE BLANK NOT FILMED 



Table C-l. Numerical Values of Parameters used for Free-Free Mode Shapes. 



Pi 




0.98250214 

1.00077731 

0.9999645 


4.73000407 


7.85320462 

10.99560783 


The mode shapes defined by the parameters in Table C-l were substituted into equations 
defining the nondimensional integrals, Eqs. (3.15) and (3.16). The numerical integration scheme 
was then used to determine the value of these integrals. 

The integrals are needed to form the constant matrices that define the dynamics of the 
flexible structure. The symbolic form of the matrices were given by Eqs. (3.33), (3.37), and 
(3.38). They are rewritten here for convenience 

[Moh = h{U) (C.5) 


[C 0 \i = 2C Va/i(m)/ 2 (m) 


(C.6) 


[K a \i = XI 2(1,1) 


(C.7) 


The numerical values for the constant mass, damping, and stiffness matrices are shown for 
the first five modes. These modes correspond to rigid body translation, rigid body rotation, and 
the first three flexible modes, respectively. The matrices are: 



(C.8) 
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0 


0 


[K 0 J=A 


500.56 

3804.53 

14,617.63 


(C.10) 


The time-varying components, due to the dynamics of the moving mass, are added to these 
matrices. The resulting matrices are used to simulate the motion of the entire system. 

C.2 THE SPATIAL DERIVATIVES OF THE SHAPE FUNCTIONS 


In the discrete formulation detailed in Section 3.2, the spatial derivatives of the shape 
functions were kept in their symbolic form. In this appendix, the required derivatives are obtained. 
Both shape functions, linear and cubic, are examined. The spatial derivatives are with respect to 
the weighting function £, which was defined by Eq. (3.85). Both the first and second spatial 
derivatives are evaluated. 


Linear Shape Function 

The linear shape function was defined by Eq. (3.87) and is rewritten here for convenience: 

Vi=(\-s)v i + sv i+1 , ru 


where V) and V /+7 are constant vectors that select the appropriate nodal values. 


The first spatial derivative, with respect to £, of Vj is 

Vi { = V M - Vi 


(C.12) 


Because Vj is linear with respect to the weighting function, the second spatial derivative is zero. 
In the formulation, an impulse force was used to artificially add the effects of the missing 
derivative. To obtain a non-zero second derivative, a cubic shape function was formed. 
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Cubic Shape Function 


The cubic shape function was defined in 


Eq. (3.92) and is rewritten here for convenience: 


V 3 =(l-3S 2 + 2S 3 )v,- + («-2^ 3 )JLv i+1 

+ (3 £ - 2 £ 3 ) Vi+ 2 + i ^ 2 + £ ) n ^‘ +3 (C.13) 


where the V% 2 and the V M vectors are constant vector that serve a similar purpose to V,- and 
V M . 

The first and second spatial derivatives of Eq. (C.13) are 

V 3{ = (-6£+ 6S 2 )Vi + (l-4£+3£ 2 )^ V i+ i 

+ (6 £ - 6 £ 2 ) V i+ 2 + (2^+3 £ 2 )l-V i+ 3 (C.14) 

V 3xt = (-6+ 12 V, + ( - 4 +6$)]Lv i+1 

+ (6 - 12 $) V i+2 + ( 2 + 6 S)J-V /+3 (C.15) 

As shown by Eq. (C.15), the second derivative of Vj is not only continuous but linear with respect 

to£ 
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APPENDIX D 

DESCRIPTION OF COMPUTER CODE 


Appendix D presents a general description of the computer code written to simulate the 
motion of the flexible-beam/moving-mass system. Five different systems were simulated: 
continuous-fixed, continuous-free, discrete-fixed, discrete-free, and discrete-multipoint-of-contact. 
Each simulation follows the same general format: the general procedure used in creating the 
simulations is explained, the differences between the systems examined are discussed and as an 
example, a listing of the code written to simulate the discrete free system is presented. 

The code was written using PRO-MATLAB (Ref. [23]), which is a product of The 
MathWorks, Inc. MATLAB is a high-performance interactive software package. It is designed for 
scientific and engineering numeric computation. In MATLAB, the problem solutions are 
expressed almost exactly as they are written mathematically. The simplicity of programming in 
PRO-MATLAB is having a matrix defined as a basic data element that does not require 
dimensioning (Ref. [23]). The simulations were performed on a UNIX-based SUN 
SPARCstation 2. 

The purpose of the code was to numerically integrate the dynamic equations describing the 
motion of the flexible-beam/moving-mass system. The parameters needed to describe the physical 
properties of the system are X, fim, and fig, all of which were described in detail in Chapter 3. 
The values of these parameters used in the simulations were listed in Chapter 4. 

A flowchart for the program is shown in Figure D-l. The required inputs are: the initial 
state vector x 0 , the initial time t 0 , the final time tf, and the nondimensional parameters. The 
desired output is the time history of the beam’s displacement. 

The constant components of the mass, damping, stiffness, and force matrices are 
determined. Because these are independent of time, they can be calculated outside of the numerical 
integration loop. After these constant matrices are known, a subroutine that performs the fourth- 
order Runge Kutta integration scheme is called, which, in turn, calls a dynamic subroutine to set 
up the matrix equation describing the system's dynamics. 

The dynamic subroutine first formulates the time-varying components of the matrices. 
These components are then added to the constant matrices that were developed earlier. The total 
mass, damping, stiffness, and force matrices are combined into the form shown in Eq. (3.32). 
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Figure D-l. Flow chart for numerical simulations. 
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The output of this dynamic subroutine is the derivative of the state vector x . This vector is 
sent back to the Runge Kutta integration routine for numerical integration. This process, Runge 
Kutta-dynamic-Runge Kutta, is repeated until the entire time history of the state vector x is known. 
Since a fourth-order Runge Kutta integration scheme is used, the dynamic subroutine is called four 
times for each time step. 

The state vector contains the time history of the modal displacements and velocities. The 
modal displacements are transformed to the natural coordinates of the beam using Eq. (3.1 1). The 
final result is the time history of the beam's displacement. 

The procedure outlined by the flowchart shown in Figure D-l represents the general format 
of the code. The following paragraphs discuss the variations used to formulate each system. 


D.l CONTINUOUS FORMULATIONS 

The two continuous formulations differ only by the mode shapes used to describe the 
beam's deformation. This difference effects the constant mass, damping, and stiffness matrices as 
these matrices are functions of the nondimensional integrals and the nondimensional frequencies, 
which themselves are a function of the mode shapes. 

For the simply-supported beam, the frequencies and the nondimensional integrals are 
determined analytically. In the code, they appear in their symbolic form. For the free-free beam, 
the nondimensional frequencies are inputed by hand according to the values dictated in Table C-l. 
The nondimenisonal integrals are evaluated using a Runge Kutta integration scheme. 

Once these matrices are developed, the two continuous formulations are almost identical 
and both follow the general procedure described above. 

D.2 DISCRETE FORMULATIONS 

The discrete formulations are a little more involved. First, the finite element mass and 
stiffness matrices are determined. Using these matrices, the mode shapes and the nondimensional 
frequencies of the system are determined. If the beam is free-free, the matrices are not altered 


129 



before the eigenvalue decomposition. For the simply-supported beam, the boundary conditions at 
the endpoints must be implemented before the eigenvalue decomposition can take place. 

These mode shapes and frequencies are next truncated to contain only the number of modes 
used in the simulation. These truncated values are used to form the constant discrete stiffness and 

damping matrices. 

To formulate the time-varying components, it is necessary to develop the shaping 
functions, which are evaluated at every time step. A different subroutine was written for the both 
the linear and the cubic shape functions. Unlike the linear shaping function, the cubic shaping 
function is also altered if there is more than one point of contact. 

The shape functions are used to formulate the time-varying components of the matrices. 
Once the matrices are known, the same procedure outlined above is followed. 

D.3 LISTING OF CODES 

A listing of the code used to solve the discrete free single point of contact system is 
presented. The cubic shape function is used to model the displacements 
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Function [q,t] = mload(n,m,us,ug,um) 

% 

% This function will simulate the motion 
% of the beam due to the moving load. 

% 

% The inputs to the system are 
% n = number of finite elements 
% m = number of modes retatined 
% us = nondimensional stiffness parameter 
% ug = nondimensional load parameter 
% um = nondimensional mass ratio 
% 

% It will first set up all the system 
% properties. 

% 

% Then it will use this property to 
% create all the variables that are 
% not dependent on time. 

% 

% It will next call integrate which will 
% numerically integrate the state vector 
% 

% In order to do this the initial state 
% vector must be set up 
% 

% After the numerical integration of 
% the state vector, the modes will be transformed 
% to the physical non-dimensional deflection 
% at the midspan of the beam. 

% 

% Set up the variables that are not dependent on time 
% 

% Compute the eigenvalues, discrete m and k matrices 

% 

[phi,w,k,mt] = pinnedf3(n); 

% 
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% Normalize the eigevector wrt to m 


% 

phin = normeig(phi,mt,n); 

% 

% Reduce to the number of modes retained 
% 

[phir.wr] = reduce3(phin,w,n,m); 

% 

% Caluculate constant compontents of K and C 
% 

[c,kr] = stiffdamp(wr,m,us); 

% 

% Now set up variables needed for integrate 
% 

% Time 
% 

tO = 0.0; 

tf = 0.010; 
delt= .001; 

% 

% Inital state vector 

% 

[lam.sig] = coef(l); 
q0 = qint(lam,sig,n); 
xO = xint(phir,q0); 
j=l+m; 

while j < 2*m + 1 
x0(j,l) = 0.0; 

j=j+i; 

end 


% Call integrate 


[tt,x] = integrate(tO,tf,delt,xO,n,m,um,ug,phir,c,kr); 
% 
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% Next do a transformation back to 
% physical coordinates 
% 

% Call trans 

% 

x=x'; 

qq = trans(phir,x,m); 
skip = 10; 

i=i; 

m=l; 

[n,j] = size(tt); 
while i < j+1 
q(:.m) = qq(:,i); 
t(m) = tt(i); 
i=i+skip; 
m=m+l; 
end 



Function [x,y,kk,mm] = pinned3f(n) 


% 

% This function will calculate 
% the mass and stiffness matrices 
% for the free free beam. 

% 

% It will also calculate the eigenvalues and eignevectors 
% of the pinned - pinned 
% 

% The program is broken up into two parts 

% 

% Part one sets up the global mass matrix 
% and the global stiffness matrix 
% 

% Part two sets up the eigenvalue problem and solves 
% for the natural frequencies of the beam. The boundary 
% conditions (if there are any) would also taken care of in this part. 

% 

% INPUTS 
% 

% n represents the number of elements to be used 
% 

% OUTPUTS 
% 

% y represents the matrix of eigenwalues that represent 
% the frequencies of the beam 
% 

% x represents the matrix of eigenvectors that will 
% will be used as the shape functions during modal 
% analysis 

% 

% But First 

% 

% Certain variables are declared 
% 

% LENGTH OF BEAM 
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1 = 1 . 0 ; 


% 

% DENSITY OF BEAM 
p=i; 

% 

% MODULUS OF ELASTICITY 
e=l; 

% 

% WIDTH OF BEAM 
w=l; 

% 

% DEPTH OF BEAM 
d=l; 

% 

% AREA MOMENT OF INERTIA OF BEAM 
im= 1; 

% 

% 

% The length of one element is 
le = 1/n; 


% PART ONE 

% 

kll = [12 6*le 
6*le 4*le A 2]; 
kl2 = [-12 6*le 
-6*le 2*le A 2]; 
k21 = [-12 -6*le 
6*le 2*le A 2]; 
k22 = [12 -6*le 
-6*le 4*le A 2]; 

% 


mil =[156 22*le 
22*le 4*le A 2]; 
ml2 = [54 -13*le 
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13*le -3*le A 2]; 
m21 = [54 13*le 
-13*le -3*le A 2]; 
m22 = [156 -22*le 
-22*le 4*le A 2]; 

% 

% Assemble the global stiffness matrix and mass matrix 

% 

k(l:2,l:2)=kll; 

k(2*n+ 1 :2*n+2,2*n+ 1 :2*n+2)=k22; 

k(2*n+l:2*n+2,2*n-l:2*n)=k21; 

k(l:2,3:4)=kl2; 

% 

m(l:2,l:2)=mll; 

m(2*n+l:2*n+2,2*n+l :2*n+2)=m22; 

m(2*n+l:2*n+2,2*n-l:2*n)=m21; 

m(l:2,3:4)=ml2; 

% 

i=2; 

while i <= ((2*n)/2) 
k(2*i- 1 :2*i,2*i-3:2*i-2)=k2 1 ; 
k(2*i- 1 :2*i,2*i- 1 :2*i)=kl l+k22; 
k(2*i- 1 :2*i,2*i+ 1 :2*i+2)=kl2; 

% 

m(2*i- 1 :2*i,2*i-3:2*i-2)=m21; 
m(2*i- 1 :2*i,2*i- 1 :2*i)=ml l+m22; 
m(2*i- 1 :2*i,2*i+l :2*i+2)=ml 2; 
i=i+l; 
end 

% 

kk=e*im*k/le A 3; 
mm = (m*le)/420; 

% 

% PART TWO 

% 

% 
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% 

% Next do an eigenvalue decomposition 
% 

[xx,y] = eig(kk,mm); 
y = diag(y); 

% 

% Next diagonalize the eigenvector matrix 
% 

plo = xx(:,2*n+l); 
p2o = xx(:,2*n+2); 

% 

alpha = l/sqrt(plo'*mm*plo); 
pin = alpha*plo; 

% 

p2n = p2o - (p2o'*mm*p 1 n)*p 1 n; 

% 

x = [xx(:,l:2*n) pin p2n]; 
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Function phin = normeig(phi,mt,n) 

% 

% This function will normalize the 
% eigenvector matrix so 
% phi'*mt*phi = I 
% 

% 

i = l; 

while i < 2*n+3 

alpha = sqrt(phi(:,i).'*mt*phi(:,i)); 
phin(:,i) = phi(:,i)./alpha; 
i=i+l; 
end 
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Function [phir,yr] = reduce3(phi,y,n,m) 

% 

% This function will reduce the 
% matrix of eigenvectors to the 
% number of modes that are 
% actually wanted in the simulation 
% versus the number of degrees of 
% freedom in the model 
% 

% The inputs are 

% phi = the normalized matrix of eigenvectors 
% y = the eigenvalues 
% n = number of elements used 
% m = number of modes wanted 
% 

% The output is 

% phir = the normalizes and reduced eigenvector matrix 
% yr = the reduced eigenvalue matrix 

% 

diff = 2*n-m+3; 

phir=phi(:,diff:2*n+2); 

yr=y(diff:2*n+2); 
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Function [c,kr] = stiffdamp(yr,m,us) 

% 

% This function will return 
% the stiffness and damping 
% matrices to be used in the 
% state vector. 

% 

% Its inputs are 
% yr = the reduced eigenvalue 
% matrix 

% m = the number of reduced 
% modes 

% us = nondimensional stiffness parameter 


j=i; 

chi = 0.01; 
while j < m + 1 
c(j j) = 2*chi*sqrt(yr(j)*us); 
kr(j j) = yr(j)*us; 

j=j+i; 

end 
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Function [tout,xout] = integrate(t<Mf,delt,xO,n,m,um,ug,phir,c,kr) 

% 

% This function will numerically integrate 
% the state vector to simulate the motion 
% of the model 
% 

% The inputs are 
% tO = inital time 
% tf = final time 
% xO = inital state vector 
% delt = delta t 
% All the rest are 

% model parameters which have 

% previosly been defined 
% 

% The output is 

% tout = the column vector containing all the times used 
% xout = the final state vector 
% 

% First initialize and set the tolerance 
% 

flag=(tf-tO)/delt; 
flag2 = 2*m; 
xout = zeros(flag,flag2); 
tout = zeros(l,flag); 

% 

% 

xout(l,l:flag2) = xO'; 
tout(l) = tO; 


xk = xO; 
t = tO; 

% 

tol = .le-06; 

% 

% Next set up the loop 
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% 

count=2; 

% 

while t < tf 

xdl = eqn3(t,xk ) n,m,um,ug,phir,c,kr); 
bl = delt*xdl; 

% 

tt = t + delt*0.5; 
xt = xk + 0.5*bl; 

xd2 = eqn3(tt,xt > n,m,um,ug,phir,c,kr); 
b2 = delt*xd2; 

% 

xt = xk + 0.5*b2; 

xd3 = eqn3(tt,xt,n,m,um,ug,phir,c,kr); 
b3 = delt*xd3; 

% 

tt = t + delt; 
xt = xk + b3; 

xd4 = eqn3(tt,xt,n,m,um,ug,phir,c,kr); 
b4 = delt*xd4; 

% 

x = xk + (l/6)*(bl + 2*b2 + 2*b3 + b4); 
% 

% Set up output vectors 
% 

t = t+delt; 

xout(count,l:flag2) = x'; 
tout(count) = t; 

% 

count = count +1; 

% 

% Get ready for next iteration 
% 

xk = x; 

% 

end 
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Function xd = eqn3(t,x,n,m,um,ug,phir,c,kr) 

% 

% This function is called by 
% integrate 
% 

% It sets up the state vector and its derivatives 

% 

% The inputs are 
% t = time 

% x = state vector at time t 

% The rest of parameters have been defined previously. 

% 

% First set up the load matrices 
% which are dependent on time 
% 

if t >=1.0 
um = 0.0; 
ug = 0.0; 
end 
% 

% Call loadf3 

% 

[v,vd,vdd] = loadf3(t,n); 

[d,e] = size(v); 
if d== 1 
v=v'; 
end 

[g,h] = size(vd); 
if g ==1 
vd=vd'; 
end 

% 

% Call loadn3 
% 

[vn,vdn,vddn] = loadn3(v,vd,vdd,phir,t); 

% 
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% 

% Next set up the mass, stiffness, damping and force matrices 
% 

% 


ma = mass(m,um,vn,t); 
da = damp(c,um,vn,vdn); 
ka = stiff(kr,um,vn,vddn); 
% fa = force(n,ug,t,phir); 
% 

% Set up the state vector 


% Without r and theta. Simply Supported Case 


%mn = ma(3:2+m,3:2+m); 

%dn = da(3:2+m,3:2+m); 

%kn = ka(3:2+m,3:2+m); 

%fn = fa(3:2+m,l); 

% 

% 

%xd(3:2+m,:) = x(5+m:2*(2+m),:); 

%xd(5+m:2*(2+m),:) = inv(mn)*(fn - kn*x(3:2+m,:) - dn*xd(3:2+m,:»; 
% 

% 

% With r and theta, Free-Free Case 
% 

xd(l:m,:) = x(l+m:2*m,:); 

% 

% With Forcing term 
% 

%xd(l+m:2*m,:) = inv(ma)*(fa - ka*x(l:m,:) - da*xd(l:m,:)); 

% 

%Without Forcing term 

xd(l+m:2*m,:) = inv(ma)*(- ka*x(l:m,:) - da*xd(l:m,:)); 
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Function [v,vd,vdd] = loadf3(t,n) 


% 

% This function will return the vectors 
% v, vd, vdd, which weights the effect of the moving 
% load on the neigboring nodes. It is a 
% function of time since the load is moving 


% The inputs are 
% t is time 

% n is number of elements used 
% for the modeling 
% 

% The outputs are 
% v cubic weighting function 
% vd first derivative of v 
% vdd second derivative of v 


% Evaluate i, and alpha 
% All of which are a function of 
% time 
% 

i = fix(t*n) + 1; 
alpha = (t*n - i + 1); 

% 

% Set up a dummy array of zeros 
% 

vl = zeros(l,2*n+2); 
v2 = zeros(l,2*n+2); 
v3 = zeros(l,2*n+2); 
v4 = zeros(l,2*n+2); 

% 

% Determine v for n less than the 
% number of elements 
% 
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if i < n+1 
vl(l,2*i-l) = 1; 
v2(l,2*i) = 1; 
v3(l,2*i+l) = 1; 
v4(l,2*i+2) = 1; 
end 
% 

% Correct last value of 
% v2 when appropriate 


if i = n+1 
vl(l,2*i-l) = 1; 
v2(l,2*i) = 1; 
end 
% 

va = (l-3*alpha A 2 + 2*alpha A 3)*vl + (alpha - 2*alpha A 2 + alpha A 3)*(l/n)*v2; 
vb = (3*alpha A 2 - 2*alpha A 3)*v3 + (-alpha A 2 + alpha A 3)*(l/n)*v4; 

% 

v = va + vb; 

% 

vc = (-6*alpha + 6*alpha A 2)*vl + (l-4*alpha+3*alpha A 2)*(l/n)*v2; 
ve = (6*alpha - 6*alpha A 2)*v3 + (-2*alpha + 3*alpha A 2)*(l/n)*v4; 

% 

vd = (vc + ve)*n; 

% 

vcc = (-6 + 12*alpha)*vl + (-4+6*alpha)*(l/n)*v2; 
vccc = (6 - 12*alpha)*v3 + (-2 + 6*alpha)*(l/n)*v4; 

% 

vdd = (vcc + vccc)*(n A 2); 

% 

v = v’; 
vd = vd'; 
vdd = vdd'; 
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Function [vn,vdn,vddn] = loadn(v,vd,vdd,phir,t) 

% 

% This function will calculate 
% vn = phir’*v 
% vdn = phir’*vd 
% vddn = phir'*vdd 

% 

% phir has already been reduced to the appropriate 
% amount of modes desired. 

% 

vn = phir’*v; 
vdn = vd'*phir; 
vddn = vdd'*phir; 
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Function ma = mass(m,um,vn,t) 

% 

% This function will calculate the mass matrix 
% 

% 

% The inputs are 
% 

% m = number of modes used 
% um = non dimensional mass parameter 
% vn = load vector 
% t = non dimensional time 
% 


% 

% 

ma = eye(m,m) + um*vn*vn'; 
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Function da = damp(c,um,vn,vdn) 

% 

% This function will calculate the 
% damping matrix 
% 

% The inputs are 

% c = diagonal structual damping matrix 
% um = non dimensional mass parameter 
% vn = load vector 
% vdn = derivative of load vector 

% 

da = c + 2*um*vn*vdn; 
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Function ka = stiff(kr,um,vn,vddn); 

% 

% This function will calculate the stiffess 
% matrix 
% 

% Its inputs are 

% kr = constant component of k 
% um = mass ratio 
% load vectors vn and vddn 
% 

fy = um*vn*vddn; 
ka = kr + fy; 
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Function x = trans(phir,xx,m) 

% 

% This function will transform the 
% modal displacements into natural 
% displacements 
% 

x = phir*xx(l:m,:); 
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